Data Science Challenge - Housing market of Santiago de Chile¶
By Andreas Kaae, s154395 & Rasmus Riber, s154060
Table of Contents
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
Introduction¶
Housing is central in the lives of almost everyone. For most people housing is responsible for a significant part of their personal economy. This notebook report focuses on the housing market in Santiago de Chile, the capitol of Chile. A dataset of housing transactions serves as the basis for the report. For each transaction, area charateristics, the type of the buying household and geographical data is avaiable.
The report consists of five main sections:
IV) Data preperation
In this section the data is described and cleaned.
V) Additional data
In this section additional data of crimes in Santiago de Chile is gathered and added to the original data. This addition is performed with the expectation of criminality having high relevance for the real estate market.
VI) Initial data analysis
In this section statistics, correlations, outliers and distributions of features are analysed.
VII) Exploratory analysis
In this section 5 research questions are analysed:
a) How is wealth related to type of household.
b) How is wealth related to other real estate features.
c) How is crime related to household characteristics.
d) Are there any geographical patterns for the real estate features?
e) How well can the real estate market be segmented based on area characteristics?
VIII) Prediction challenge
This section will adress the research question of predicting the household type that most likely will buy a specific real estate, by appliying machine learning tools.
IX) Conclusion
This section will conclude the report by stating the key findings.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
df = pd.read_excel('Households_Data.xlsx')
df_coords = df[['DirCoordX','DirCoordY']] #Saved for later use
df.head()
df.shape
The data set describes details about households in Santiago de Chile with details about household transactions as well as the charactaristics about the household which bought the house. The main author of the data is Tomas Cox Oettinger, a PhD student from Chile. The dataset has 18264 observations and 30 features.
| Feature | Description | Type |
|---|---|---|
| Hogar | ID of the household who bought the real estate | Categorial |
| EXCLUIR | Whether or not the observation should be in the dataset | Binary |
| paraValidacion | Sample for validaiton | Binary |
| Sector | Exogenours sector | Categorical |
| DirCoordX | UTM coordinat X | Continuous |
| DirCoordY | UTM coordinat Y | Continuous |
| MontoArr | Monthly dwelling rent (Chilean Peso) | Continuous |
| IngresoHogar | Household income (Chilean Peso) | Continuous |
| Factor | Expansion factor of household survey | Continuous |
| AnosEstudio | Years of formal education of head of household | Integer |
| CLASE | 1: less than 12 years, 2: between 12 and 15 yeard, 3: more than 15 years | Categorical |
| Sup_Prom_Constr_Hab_EOD | Average built surface of residential units in zone | Continuous |
| Calid_EOD_norm_inv | Average Construction Quality of residential units in zone | Continuous |
| DensConstr_EOD | Built density in zone (built suerface divided by zone surface) | Continuous |
| Dist_est_MEtro_MC_KM | Distance to nearest subway station (km) | Continuous |
| Dist_salida_Autop_MC_KM | Distance to nearest highway exit (km) | Continuous |
| Tiempo_Com_Stgo | Travel time to CBD (min) | Continuous |
| Ingreso_Promedio_Zone_MM | Average zonal income (Million Chilean Peso) | Continuous |
| Acc_Comercio_tpte_pub | Accesibility to Commerce in transit | Continuous |
| Acc_Educacion_tpte_pub | Accesibility to Educaition in transit | Continuous |
| Acc_Habitacion_tpte_pub | Accesibility to Habitational in transit | Continuous |
| Acc_Industria_tpte_pub | Accesibility to Industry in transit | Continuous |
| Acc_Servicios_tpte_pub | Accesibility to Services in transit | Continuous |
| Acc_Comercio_tpte_auto | Accesibility to Commerce in car | Continuous |
| Acc_Educacion_tpte_auto | Accesibility to Educaition in car | Continuous |
| Acc_Habitacion_tpte_auto | Accesibility to Habitational in car | Continuous |
| Acc_Industria_tpte_auto | Accesibility to Industry in car | Continuous |
| Acc_Servicios_tpte_auto | Accesibility to Services in car | Continuous |
| CLUSTER7 | Cluster clasification (7 clusters) | Continuous |
| CLUSTER2 | Cluster clasification (2 clusters) | Continuous |
Coordinate transformation
The features 'DirCoordX' and 'DirCoordY' are listed in UTM format. To locate the house holds in Santiago de Chile, the coordinates are converted to WGS84 format. WGS84 is the standard format applied for maps and is applied by google maps as well. The package utm is applied to convert UTM coordinates to the format of the standard coordinates. "19, E'" indicates the UTM area of Santiago de Chile.
import utm
df['Latitude'], df['Longitude'] = utm.to_latlon(df['DirCoordX'].values, df['DirCoordY'].values, 19, 'E')
df.loc[:,['DirCoordX','Latitude','DirCoordY','Longitude']].head()
During the conversion the orginial coordinate 'DirCoordY' is apparently changed, hence the orginial coordinate columns are reloaded.
df['DirCoordX'] = df_coords['DirCoordX']
df['DirCoordY'] = df_coords['DirCoordY']
df.loc[:,['DirCoordX','Latitude','DirCoordY','Longitude']].head()
The new columns 'Latitude' and 'Longitude' (in WGS84 coordinates) are now a part of the dataset 'DirCoordX' and 'DirCoordY'.
Removal of features and missing values
In this section redundant features and missing-values will be removed. Not all features are relevant for the challenge and are therefore removed for the sake of simplicity. Initially the following columns are excluded from the data set:
- EXCLUIR - is a column specified by the author stating whether he excluded the row or not. As no reason for the exclusions is providede the column is removed.
- Hogar - is the household ID and is therefore set as index.
- paraValidaction - is a column generated by the author of the data, hence it is removed.
- Factor - is the weight of the kind of familiy in the population and deemed irrelevant for the challenge.
- CLUSTER7 and CLUSTER2 - are cluster clasification columns which are made by the data author. The clusters are based on the other features and are therefore removed.
- AnosEstudio - This is the feature which CLASE is generated from and is therefore removed.
# Delete columns
exclude_columns = ['EXCLUIR','paraValidacion','Factor','CLUSTER7','CLUSTER2','AnosEstudio'] # Columns to be excluded
df = df.drop(exclude_columns,1)
# Set index
df = df.set_index(['Hogar'])
In the feature CLASE there are 4 different categories: '-99', '1', '2' and '3'. The category '-99' means that the respondant of the survey did not answer the question and these observations are therefore excluded.
df = df[df['CLASE'] != -99]
The two categorical features CLASE and Sector are set as categorical variables.
df['CLASE'] = df['CLASE'].astype('category')
df['Sector'] = df['Sector'].astype('category')
Missing values in the dataset can cause problems, hence the data is checked for 'NaN' or similar entries. The isnull() does excatly this, if there is a 'NaN' it will return 'True'. Thus summing over the whole dataframe detects any missing values.
df.isnull().sum().sum()
As the result is 0, there are no NaN-values in the dataset.
The dataset have now been described and irrelavant entries have been removed thus reducing the dataframe to 17876 observations and 25 features.
df.shape
In this section additional data of crimes in Santiago de Chile is gathered and added to the transaction data. This addition is performed with the expectation of criminality having high relevance for the real estate market.
Data source - https://en.datachile.io/about/data
The data is avaiable in a package of three sets. Crime records (felonies.csv), geographic data (comunas_id.csv) and descriptions of crime types (categories_id.csv). In this section the three files are merged.
The crime data is loaded.
crimes = pd.read_csv('felonies.csv')
crimes.year = pd.to_datetime(crimes.year, format='%Y').dt.year
print('Start year:',min(crimes.year),' End year:',max(crimes.year))
print('Number of communes:',len(crimes.comuna_datachile_id.unique()))
crimes.head()
The felony_id.csv data consists of 4 features, and contains observations from 2005 to 2016.
| Attribute | Description | Type |
|---|---|---|
| year | The year the crime was conducted | Date |
| comuna_datachile_id | The commune id (similar to a zip) | Categorical |
| felony_id | Type of crime | Categorical |
| cases | The number of crimes | Integer |
The commune data is loaded.
communes = pd.read_csv('comunas_id.csv')
communes = communes.drop(['comuna_customs_id','comuna_tax_office_id'], axis = 1)
print(communes.shape)
communes.head()
The comunas_id.csv data consists of 4 features for 348 communes.
| Attribute | Description | Type |
|---|---|---|
| region_id | ID of the region | Categorical |
| region_name | Name of the region | Categorical |
| comuna_datachile_id | ID of the commune | Categorical |
| comuna_name | Name of the commune | Categorical |
Initially it was decided to look at the Santiago de Chile metropolitan region only by:
It reduced the amount of communes to 54. But it turned out that some real estate transactions were recorded outside this region. Hence the entire data is kept for now.
The communes data is sorted by the comuna_datachile_id to do a quick inspection. Afterwards communes is merged withe the crime data.
#communes = communes[communes['region_name'] == 'Región Metropolitana Santiago']
print(communes.shape)
communes = communes.sort_values(by='comuna_datachile_id')
communes.head()
crimes = crimes.merge(communes,how='right',on='comuna_datachile_id')
print(crimes.shape)
crimes.head()
Merging the data results in the crime dataframe including region id and name as well as commune name. There is crime data for 346 out of 348 communes.
The data is checked for missing values and erros.
crimes.isnull().sum()
crimes[crimes.felony_id.isnull()]
Since only 2 observation are lacking information, these are discarded.
crimes = crimes[np.isfinite(crimes.felony_id)]
crimes['felony_id'] = crimes['felony_id'].astype(int)
print(crimes.shape)
crimes.head()
The data set describing the crime types is loaded. It consists of the 14 felony ids and the felony desriptions in Spanish.
categories = pd.read_csv('categories_id.csv')
categories.head(14)
The dataframe is converted into a dictionary to easily insert the descriptions in the crimes dataset.
categories = categories.set_index('felony_id')
categories = categories.to_dict('dict')
crimes['felony_es'] = crimes['felony_id'].map(categories['felony_es'])
crimes.head()
crimes.isnull().sum()
The crime data sets are now combined.
In this section the crime data is connected to the dataframe describing the households. This is done in three steps by using the commune id and the corresponding geographic location.
1) The first step it to load the geographic data (ZonasEOD2012.dbf) for the communes in Chile.
2) The second step identifies the commune id of the transaction real estates.
3) The third step adds the criminality data to the transaction data via the commune id.
1) Load geographic data for the communes in Chile¶
Initially the geographic data is loaded. To work with geographic data the packages geopandas and shapely.geometry are used.
import geopandas as gdp
from shapely.geometry import Point, Polygon
#Load the data provided about the zones:
zones = gdp.read_file('Zones/ZonasEOD2012.dbf')
# Define crs zone
crs_zones = zones.crs
zones['comuna_datachile_id'] = zones['ID'].astype(int)
zones = zones[['comuna_datachile_id','geometry']]
zones.head()
print(type(zones))
print(crs_zones)
The loaded dataframe is a geopandas dataframe and the code used for the coordinates is epsg:32719. This 'key' is used below to specify which coordinate type is used.
Secondly the commune data is added to the real estate transcation data via the coordinates. For each real estate coordinate, a geographic point is zipped together from its x (DirCoordX) and y (DirCoordY) coordinates.
df = gdp.GeoDataFrame(df, crs = crs_zones, geometry = [Point(x, y) for x, y in zip(df.DirCoordX, df.DirCoordY)])
df.head()
fig,ax = plt.subplots(figsize=(15,15))
# Plot the communas
zones.geometry.plot(ax=ax, alpha = 1, color = 'grey')
# Plot the household observations
df.geometry.plot(ax = ax, markersize = 1, color = 'tab:blue')
plt.show()
The map illustrates the communes (grey areas) and the househould observations (blue dots). Finally, the comunes in the transaction data and the crime data is compared to check how well they match.
print('Communes in crime data:',len(crimes.comuna_datachile_id.unique()))
print('Communes in zone data:',len(zones.comuna_datachile_id.unique()))
print('Comumunes in both crime and zone data:',len( set(crimes.comuna_datachile_id.unique()) & set(zones.comuna_datachile_id.unique()) ))
2) Identifyng commune id of the transaction real estates¶
In this section the commmune is added to the transaction data. The code is outcommented due to long run time. The code is kept to show how the data was prepared and saved as a csv. Eventualy the csv is loaded to save run time.
The section is inspired by this page https://stackoverflow.com/questions/48097742/geopandas-point-in-polygon
Each hogar is assigned to a commune by matching its coordinates to the zones data as illustrated below:
#Lists for polys
#c_ids = zones.comuna_datachile_id.tolist()
#geos = zones.geometry.tolist()
#Polys
#polys = gdp.GeoSeries({})
#for row in range(len(c_ids)):
# polys[str(c_ids[row])] = geos[row]
#points
#_pnts = df.geometry.tolist()
#pnts = gdp.GeoDataFrame(geometry=_pnts, index=df.index)
#pnts.head()
#Assign points to comuna - OBS OBS OBS takes 1-2 minutes
#pnts = pnts.assign(**{key: pnts.within(geom) for key, geom in polys.items()})
#pnts.head()
#print('Total transactions:',len(df))
#print('Transactions assigned to 0 communes:',len(pnts[pnts.select_dtypes(include=['bool']).sum(axis=1)==0]))
#print('Transactions assigned to 1 commune',len(pnts[pnts.select_dtypes(include=['bool']).sum(axis=1)==1]))
#print('Transactions assigned to 2 communes',len(pnts[pnts.select_dtypes(include=['bool']).sum(axis=1)==2]))
#print('Transactions assigned to 3+ communes',len(pnts[pnts.select_dtypes(include=['bool']).sum(axis=1)>2]))
The 43 transactions that aren't assigned to a commune is left out. The 20 transaction assigned to 2 communes are assigned to the last. A new column indicationg the comunas_id is addded.
#pnts['comunas_id'] = np.nan
#pnts.head()
The column is filled as per below.
#row_failures = 0
#for row in range(len(pnts.index)):
# try:
# pnts.comunas_id.iloc[row] = pnts.columns[(pnts == True).iloc[row]][0]
# except:
# row_failures += 1
# print('Row',row,' has no Trues.')
#row_failures # =43
#export_csv = pnts.to_csv('pnts.csv', index = True, header=True)
Here the table was saved. Going ahead the table is loaded as per below:
pnts = pd.read_csv('pnts.csv', usecols=['Hogar','comunas_id'])
pnts['comunas_id'] = pnts['comunas_id'].apply(lambda x: int(x) if x == x else "") #deals with float to int and deal with NaN
pnts.head()
3) Merging criminality and transaction data by commune id¶
Finally the comunas_id are merged to the transaction data set applying Hogar as the key. The index is reset and will later be set to Hogar again.
# Reset the index
df=df.reset_index()
df.head()
#Merge on Hogar to get the comunas_id column.
df = df.merge(pnts,how='left',on='Hogar')
df.head()
Finally the amount of crime (sum of all cases) in the commune is merged to the real estate transaction date. The initial intend was to categorise crime by severity but since the data turned out to be missing for the city center, hence this level of detail has been left out.
The crime data is grouped by commune and merged on commune ID to the transaction data, as illustrated below:
# Change heading of column
crimes = crimes.rename(columns={"comuna_datachile_id": "comunas_id"})
crimes.head()
# Group by communa ID
crimes_df = pd.DataFrame(crimes.groupby('comunas_id').cases.sum())
crimes_df = crimes_df.reset_index()
crimes_df.head()
The two dataframes are now meged so the household transactions dataframe gets the amount of cases within the communas_id.
df = df.merge(crimes_df,how='left',on='comunas_id') #Merge on Hogar to get the comunas_id column.
df.head()
print('Total communes:',len(df.comunas_id.unique()))
print('Communes with crime records:',len(df[df.cases > 0].comunas_id.unique()))
print('Total transactions:',len(df))
print('Real estates with crime records:',len(df[df.cases > 0]))
One third of the transactions are in communes where crime data is avaiable. For the communes without crime date, the mean amount of crime in all communes is applied. This will allow modelling and give a neutral effect in models.
At the same time DirCoordX, DirCoordY and geometry are removed as they are explaiend by Longitude, Latitude and communas_ID and are no longer needed.
df['cases'] = df['cases'].fillna(df['cases'].mean())
coords_and_geom = df[['Hogar','DirCoordX', 'DirCoordY','geometry']] #Save for potential use later on.
df = df.set_index('Hogar')
df = df.drop(['DirCoordX', 'DirCoordY','geometry'], axis=1)
df.head()
A more throughout approach of filling out the missing data, could be to take a mean of the neighbouring communes only. However, this level of detail is not deemed to be key for this project.
Another approach to processing the crime data could be to keep the crimes seperated by type or group them by severity. Again this level of detail is not key for this project. The intend of adding the crime data is purely to get insight into any patterns between the level of crime and the type of households in an area.
Finally the features cases is renamed to crimes
df = df.rename(columns={"cases": "crimes"})
Initial data analysis¶
In this section statistics, correlations, outliers and distributions of features are analysed.
Statistical overview¶
A quick statistical overview of the features is listed below.
df.describe()
As most attributes will be investigated further, only statistics for a single attribute is explained here to illustrate the understanding.
The mean monthly rent 'MontoArr' is 186k Chilean Pesos and the IQR is [110k, 200k] Chilean Pesos. The minium and maximum monthly rents are 10k and 500k Chilean Pesos respectively. Below a boxplot and the distribution for 'MontoArr' is shown.
## Set subplot
fig, axs = plt.subplots(1, 2, figsize=(15,5))
## Plots
axs[0].boxplot(df.MontoArr);
axs[1].hist(df.MontoArr,bins=100);
axs[0].set_title("Boxplot of the feature 'MontoArr'")
axs[1].set_title("Distribution of the feature MontoArr'")
fig.tight_layout()
plt.show()
print('Households with MontoArr above 1.5M:',len(df.MontoArr[df.MontoArr>1500000]))
The median monthly dwelling is 150k Chilean Pesos. The median is lower than the mean indicating that a minority of the households pays a much higher rent than the majority. This is also illustrated by the many outliers in the boxplot and by the historgram having a positive skew.
21 households have significantly high monthly rents (1.5M). To get a better overview the MontoArr data is illustrated again without these.
## Set subplot
fig, axs = plt.subplots(1, 2, figsize=(15,5))
## Plots
axs[0].boxplot(df.MontoArr[df.MontoArr < 1500000]);
axs[1].hist(df.MontoArr[df.MontoArr < 1500000],bins=100);
axs[0].set_title("Boxplot of the feature 'MontoArr'")
axs[1].set_title("Distribution of the feature 'MontoArr'")
fig.tight_layout()
plt.show()
The majoroty of the households, the ones in the IQR, all pay a very similar rent. A significant part of the households pays a much higher rent.
Correlation check
In this section the correlations between features are looked into to determine if any should be removed. Since Clase is a categorical feature it is not considered in this correlation check. Instead correlation between the continous features are investigated.
It is desirable that explanatory features are uncorrelated, since highly correlated features will not add any relevant information to exaplin the target, but instead increase the dimensionality of the dataset and the computational time when analysing. It is however desirable that the explanatory features are correlated with the target feature. Since the target is categorical it does mathematical not make sense to look at the correlation between it and features, however this is something practiced when doing data analysis, the same is therefor done here. CLASE is therefor set back to being integer of type.
To get a better understanding of the correlations, a Person correlation matrix is plotted which shows the correlation between the continous features.
- A value close to 0 means there is little correlation (0 meaning excatly no correlation)
- Values close to 1 menas there is a strong positiv correlation
- Values close to -1 means there is a strong negative correlation
# Reset CLASE as continous
df['CLASE'] = df['CLASE'].astype('int64')
#Using Pearson Correlation
plt.figure(figsize=(20,10))
cor = df.corr()
sns.heatmap(cor, annot=True, cmap=plt.cm.Reds)
plt.show()
A key observation from the correlation matrix is that there is a large area with dark red colours from the feature 'Acc_Comercio_tpte_pub' to the feature 'Acc_Servicios_auto'. These features are highly positive correlated and will be evaluated further. Furthermore 'Dist_est_Metro_MC_KM' is highly negativ correlated with 'Acc_Educacion_auto' and 'Acc_Habitacion_auto' and should also be investigated.
All the features starting with 'Acc' and all the ones which ends with 'pub' has a correlation higher than 0.8 and most have higher than 0.95 with each other. The same is true for the features which starts with 'Acc' and ends on 'auto'. Therefor all but one feature can be removed within each of these two groups, without loosing significant explanatory power. As 'Acc_Servicios_tpte_pub' and 'Acc_Servicios_auto' are the two features which are highly correlated with the target variable these two features are selected to keep. By removing these the negativ correlations for 'Dist_est_Metro_MC_KM' are also dealt with.
The original dataframe is saved as it will be used later and the columns specified are removed.
# Set CLASE as categorical
df['CLASE'] = df['CLASE'].astype('category')
# Save the orginal cleaned data set
df_all = df
# Columns to exclude
exclude_columns = ['Acc_Comercio_tpte_pub','Acc_Educacion_tpte_pub', 'Acc_Habitacion_tpte_pub','Acc_Industria_tpte_pub',
'Acc_Comercio_auto','Acc_Educacion_auto', 'Acc_Habitacion_auto', 'Acc_Industria_auto','Acc_Servicios_auto']
df = df.drop(exclude_columns,1)
In this section the data is being investigated for outliers. Outliers with possible values will be kept in the data, but outliers with impossible data will be removed. Intially the distribution of each attribute in the data is displayed with a boxplot to identify outliers. Here the features communas_id, Sector and CLASE are left out as these are categorical featurs.
fig, axs = plt.subplots(4, 4,figsize=(15,10))
i = 0
j = 0
cont_feats = ['MontoArr', 'IngresoHogar', 'Sup_Prom_Constr_Hab_EOD', 'Calid_EOD_norm_inv',
'DensConstr_EOD','Dist_est_Metro_MC_KM', 'Dist_salida_Autop_MC_KM',
'Tiempo_Com_Stgo','Ingreso_Promedio_Zona_MM', 'Acc_Servicios_tpte_pub',
'Latitude','Longitude', 'crimes']
for att in range(len(df[cont_feats].columns)):
axs[i, j].boxplot(df[cont_feats][df[cont_feats].columns[att]])
axs[i, j].set_title('%s' %df[cont_feats].columns[att])
if j == 3:
i+=1
j=0
else:
j+=1
plt.tight_layout()
In the boxplot the upper and lower whiskers are defined as (q75+1.5IQR) and (q25-1.5IQR) respectively. Per definition, an observation outside the whisker range is considered an outlier.
Several features in the dataset have outliers. The outlier statistics for each feature is calculated below.
for att in df[cont_feats].columns:
Q1 = df[att].quantile(0.25)
Q3 = df[att].quantile(0.75)
IQR = Q3 - Q1
diff = df.shape[0] - df[(Q1-1.5*IQR <= df[att]) & (df[att] <= Q3+1.5*IQR)].shape[0]
under = len(df[Q1-1.5*IQR >= df[att]])
over = len(df[Q3+1.5*IQR <= df[att]])
if diff != 0:
print('-------',att,'-------')
print('Total outliers:',diff)
print('Total outliers %:', round(diff/df.shape[0]*100,1) )
print('Outliers above:',over)
print('Outliers below:',under)
else:
print('-------',att,'-------')
print('No outliers')
Selected attributes with a potential outlier issue are discussed here.
Wealth
The 'MontoArr' distribution clearly indicates that the upperclass is a minoirty since 9.4% of households have monthly dwellings as high outliers. Particularly two monthly dwelling are way beyond the rest. It is decided to remove these observations as they are twice as high as the second highest and do not describe the housing market very well, since few real estates will have a rent of this magnitude.
The 'IngresoHogar' distribution clearly indicates that the upperclass is a minoirty since 5.9% of households have income as high outliers. No incomes are far beyond or below all others.
The 'Ingreso_Promedio_Zona_MM' distribution is similar to that of 'MontoArr' and 'IngressoHogar' with 10.6% high outliers indicating that 10.6% of the zones in the city has a average income as high outlier.
Construction Density
The 'Sup_Prom_Constr_Hab_EOD' dsitribution has a compact IQR around 50 and 14.1% outliers with 3/5 being high outliers. This indicates that the average built surface by residential units is quit stable in 85.9% of the city. None of the outliers are far beyond or below all others.
The 'DensConstr_EOD' distribution has a comapct IQR around 0.5 and 13.6% high outliers. The attribute indicates built density in zone (built surface divided by zone surface), hence values above 1 indicates buildings with multiple floors.
Transport
Both 'Dist_est_Metro_MC_KM' and 'Dist_salida_Autop_MC_KM' have compact IQRs and no outliers far beyond the rest.
The 'Tiempo_Com_Stgo' attribute indicates the travel in minutes to CBD. Four observations have negative travel times which are impossible values. These observations are excluded from the data.
Criminality
A single commune had 500k crime cases from 2005 to 2016. This observation is far beyond the second highest and translates into a daily average of 125 crimes. The commune with the second highest level had 300k crimes from 2005 to 2016 translating into a daily average of 75 crimes. A similar approach as used for removing outliers in 'MontoArr' could have been used, but in this case the amount of crimes are for a whole area and do not differ for each observation. No observations are therefor removed because of the amount of crimes. The high level of crime is doubtful, however this can be possible with 5.6 million in habitants, and the source is deemed valid, hence this is not excluded.
#Exclude the highest 'MontoArt'
print(df.shape[0])
df = df[df['MontoArr']<4000000]
print(df.shape[0])
#Exclude Tiempo_Com_Stgo < 0
df = df[df.Tiempo_Com_Stgo > 0]
df.shape[0]
In total 6 observations are removed thus reducing the number of observations to 17870.
In this section the distribution of the target feature is checked.
## Define plot size and style
plt.rcParams['figure.figsize'] = (5, 3)
plt.style.use('ggplot')
## Bar plot
plt.bar(['CLASE = 1','CLASE = 2','CLASE = 3'], df.groupby(['CLASE']).count()['MontoArr'].values)
plt.show()
There are fewer observation of CLASE = 3 than the rest. This is something that has to be considered when training and testing a model to do prediction. It will however also depend on how the observations are distributed in the test and training sets. But as the test/train split is decided in the assignment to be tested on the city center and trained on the surrounding observations it is simply something to be aware of.
Exploratory analysis¶
This section will explore the data, to reveal valuable insights for real estate companies in Santiago de Chile.
a) Wealth and household type¶
This section explores how wealth is related to the type of household. Below the IngresoHogar (Household income) and MontoArr (monthely rent) are plotted against the 'CLASE'.
## Set subplot
fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(8,3))
## Scatterplots
ax1.scatter(df.CLASE,df.IngresoHogar)
ax1.set_title("IngresoHogar")
ax2.scatter(df.CLASE,df.MontoArr)
ax2.set_title("MontoArr")
fig.tight_layout()
plt.show()
It can be seen that the observations with highest income are the observations which are in the CLASE = 3, and thereby also the ones with the longest educations. It can further be seen that there are more observations with high rents in CLASE = 3. To see more clearly how this is distributed between the different CLASE's histograms are plotted. To better visualize the distributions the x-axis is limited, so some of the observations are not within the plot.
plt.rcParams['figure.figsize'] = (15, 5)
plt.style.use('ggplot')
fea = ['IngresoHogar','MontoArr']
limit = [5*10**6, 1250000] # limits for x-axis
for i in range(0,len(fea)):
x = df[df['CLASE'] ==1][fea[i]]
y = df[df['CLASE'] ==2][fea[i]]
z = df[df['CLASE'] ==3][fea[i]]
plt.hist([x, y, z], bins=100, label=['CLASE = 1', 'CLASE = 2','CLASE = 3'])
plt.legend(loc='upper right')
plt.title(fea[i])
plt.xlabel("Chilean pesos")
plt.ylabel("Amount of observations")
plt.xlim(0,limit[i])
plt.show()
plt.show()
The two bar plots back up the patterns showed in the scatter plots. The CLASE = 1 has most observations when the household income and monthly rent are lowest. As the income and rent gets bigger (moves towards the right in the plots) it becomes the CLASE = 2 which has most observations and for the highest incomes and rents CLASE = 3 has the most observations. It thus seems that the longer an education, the higher the income and the more expensive the rent will be.
b) Wealth and area characteristics¶
This section looks into how the wealth is related to area characteristics to see how wealth determines what kind of real estate one buys. Below scatter plot for selected features are plotted.
fig, axes = plt.subplots(2, 3, figsize=(15,8))
sns.scatterplot(df['IngresoHogar'], df['Dist_est_Metro_MC_KM'], data=df, ax = axes[0,0])
sns.scatterplot(df['IngresoHogar'], df['Acc_Servicios_tpte_pub'], data=df, ax = axes[0,1])
sns.scatterplot(df['IngresoHogar'], df['Tiempo_Com_Stgo'], data=df, ax = axes[0,2])
sns.scatterplot(df['IngresoHogar'], df['crimes'], data=df, ax = axes[1,0])
sns.scatterplot(df['IngresoHogar'], df['Dist_salida_Autop_MC_KM'], data=df, ax = axes[1,1])
sns.scatterplot(df['IngresoHogar'], df['Sup_Prom_Constr_Hab_EOD'], data=df, ax = axes[1,2])
fig.tight_layout()
plt.show()
The plots show the following tendencies:
- Dist_est_Metro_MC_KM: Many people with low income have long distances to the nearest metro station. No high income people live more than 10 km away from a metro station.
- Acc_Servicios_tpte_pub: There is no clear patterns for income and access to services with public transport.
- Tiempo_Com_Stgo: Similar to the insights from the metro stations, more people with low income have longer travel time to the city center. No high income people have long travel time to city center.
- crimes: The crime levels are low for the high income citizen, for some low income households the crime rate is high.
- Dist_salide_Autop_MC_KM: Many low income households have longer distances to the highway, and no high income households have long distances to the highway.
- Sup_Prom_Constr_Hab_EOD: There is no clear relationship between the average built surface of residential units in the zone and the wealth.
From this analysis it is concluded that wealth has a role in determining the characteristics of what kind of real estate a household will buy.
c) Crime¶
This section explores whether the level of crime in a commune has an impact on the type of households in the commune.
First step is to investigate whether the level of crime is different accros coomunes. The first hypothesis is that there is significant difference between the level of criminality between communes.
c_id = [str(i) for i in df.groupby('comunas_id').mean().sort_values(['crimes'],ascending=False).index]
crim = df.groupby('comunas_id').mean().sort_values(['crimes'],ascending=False).crimes
plt.bar(c_id,crim)
plt.xlabel('Comuna_id')
#plt.xticks(np.arange(len(c_id)), c_id, rotation='vertical')
plt.ylabel('Crimes from 2005 to 2016')
plt.show()
There is clearly some coomunes with much higher levels of crimininality than others. The main plateu indicates the communes reciving the mean value in the data preperation. The communes with the highest level of criminiality have much higher levels than the mean and the least criminal communes.
For the rest of this section, communes where missing values for crime have been replaced by the mean amount of crime, are excluded, to get the true insight for the real estate agents.
Checking the distribution without the mean replacemented communes, the bar chart below, indicates a steep decline in crimes from the most criminal communes to the least criminal communes. Since, the observations with NaN's for crimes were replaced by the mean, a few more observations has been removed, hence the mean +-1 is catching this skew.
crim_mean = df.crimes.mean()
df_c = df[(df.crimes < crim_mean-2) | (df.crimes > crim_mean+2)]
#Since the Nans were replaced by the mean a few more observations has been removed, hence the mean +-1 is cathing this skew.
c_id = [str(i) for i in df_c.groupby('comunas_id').mean().sort_values(['crimes'],ascending=False).index]
crim = df_c.groupby('comunas_id').mean().sort_values(['crimes'],ascending=False).crimes
plt.bar(c_id,crim)
plt.xlabel('Comuna_id')
#plt.xticks(np.arange(len(c_id)), c_id, rotation='vertical')
plt.ylabel('Crimes from 2005 to 2016')
plt.title('The criminal level for communes with original crime data')
plt.show()
k = 50
plt.rcParams['figure.figsize'] = (15,5)
plt.style.use('ggplot') #or 'default'
plt.bar(c_id[:k],crim[:k])
plt.xlabel('Comuna_id')
plt.xticks(np.arange(len(c_id[:k])), c_id[:k], rotation='vertical')
plt.ylabel('Crimes from 2005 to 2016')
plt.title('The %i most criminal communes' %k)
plt.show()
plt.bar(c_id[-k:],crim[-k:])
plt.xlabel('Comuna_id')
plt.xticks(np.arange(len(c_id[-k:])), c_id[-k:], rotation='vertical')
plt.ylabel('Crimes from 2005 to 2016')
plt.title('The %i least criminal communes' %k)
plt.show()
The 50 most criminal communes experience far more crime than the 50 least criminal communes. The hypthesis is confirmed, theres is significant difference between the level of criminality between communes. With such a significant difference in crime between communes it might impact the types of households living there, hence the exploration continues.
Second step is to search for patterns between the ratio of household types and the level of criminality in the communes. Intially the hypthesis is:
In the top 50 most criminal communes, the ratio of 'low' and 'medium' househould types (CLASE - 1,2) is higher, and the ratio of 'top' household type (CLASE - 3) is lower, than on the 50 least criminal communes.
df_topk= df_c[df_c.comunas_id.isin(list(map(int, c_id[:k])))]
df_topk= df_topk[['CLASE','comunas_id']]
g1 = df_topk.groupby(['CLASE']).count()
df_botk = df[df.comunas_id.isin(list(map(int, c_id[-k:])))]
df_botk = df_botk[['CLASE','comunas_id']]
g2 = df_botk.groupby(['CLASE']).count()
print('Amount of house holds in top 50 criminal communes:',len(df_topk))
print('Amount of house holds in the 50 least criminal communes:',len(df_botk))
# Link to color names - https://matplotlib.org/examples/color/named_colors.html
fig, axes = plt.subplots(1, 3, figsize=(20,6))
fig.suptitle('The household type split...', fontsize=16)
axes[0].set_title('...in the 50 most criminal communes')
g1.plot.pie(y='comunas_id',ax=axes[0],autopct='%1.0f%%',colors = ['aliceblue', 'deepskyblue','royalblue'])
axes[0].set_ylabel('')
axes[0].legend(['Clase 1','Clase 2','Clase 3'])
axes[1].set_title('... in the 50 least criminal communes')
g2.plot.pie(y='comunas_id',ax=axes[1],autopct='%1.0f%%',colors = ['aliceblue', 'deepskyblue','royalblue'])
axes[1].set_ylabel('')
axes[1].legend(['Clase 1','Clase 2','Clase 3'])
axes[2].bar(np.arange(len(['1','2','3'])), [-10,-5,30], width=0.9, color=('aliceblue', 'deepskyblue','royalblue'))
plt.sca(axes[2])
plt.xticks(np.arange(len(['1','2','3'])), ['1', '2', '3'])
plt.ylabel('Difference in CLASE ratio [%]')
plt.xlabel('CLASE')
plt.title('Difference')
plt.show()
Theres is approximatly 1000 observations in top 50 and bottom 50, hence the sample size is deemed valid for the this test.
The ratio of 'lowest' and 'medium' household types are 10% and 5% less in the least criminal communes compared to the most criminal, respectively. The ratio of the 'top' househould type is 30% higher in the in the least criminal communes compared to the most criminal. The hypthesis in hereby confirmed. With the confirmation of the hypothesis, an even more throughout analysis is worth conducting, hence the exploration continues.
The third step is to check whether the difference in type of household is valid throughout the communes.
The hypothesis is that there is a correlation between the average household level and the level of criminality in a commune.
The average househould type is here defined as the mean level e.g. two CLASE 1, one CLASE 2 and one CLASE 3 would yield (1+1+2+3)/4 = 1.75.
df_c = df_c.astype('int')
mean_type = df_c.groupby('comunas_id')['CLASE','crimes'].mean()
mean_type = mean_type.sort_values('crimes',ascending=False)
mean_type['x'] = np.arange(len(mean_type))
mean_type = mean_type.reset_index()
ax = mean_type.plot(x='x',y='CLASE', kind='scatter')
ax.set_xlabel("Criminality ranking of communes (1 is most criminal)")
ax.set_ylabel("Average household");
ax.axhline(y=mean_type.CLASE.mean());
The blue dots indicate the average household for each commune, ranking from the most criminal commune to the left and the least criminal to the right. The red line indicates the average household level across all commune averages.
It is clearly evident that the type of household has no correlation with criminal ranking of the commune the household lives. The hypothesis is rejected.
The communes are investigated further to check the acutal relationship between average househould type and actual level of crime. The hypothsesis is that there is no linear relationship between the actual level of crime and the average household type in communes.
ax = mean_type.plot(x='crimes',y='CLASE', kind='scatter')
ax.set_xlabel("Crimes and household type pr. commune")
ax.set_ylabel("Average household");
There is no linear relationship between the actual level of crime and the average household type in a commune. Hence, the hypothesis is confirmed. Unfortunatly this implies that the crime data is of low value for the real estate agency.
This section has revealed some valuable insights for the real estate company. Apparently, criminal communes does not attract less household of high CLASE.
d) Geography¶
This section explores geographical patterns for the real estate market. The main hypothesis is that the type of household is divided geographificaly. If this it the case the household features and area characteristics are explored as well.
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig, ax = plt.subplots(figsize=(20,9))
cm = plt.cm.get_cmap('RdYlBu_r')
sc = ax.scatter(df.Longitude,df.Latitude, s=1, c=df['CLASE'], cmap=cm)
ax.set_title('CLASE by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
cbar = fig.colorbar(sc,ax=ax,ticks=[1, 2, 3])
cbar.ax.set_yticklabels([1, 2, 3]);
The map clearly illustrates a significant division of household type by geography, hence the hypothesis is confirmed. Northeast is completely dominated by household type 3. Mid-east is mixed. Southeast is dominated by type 2. Southwest is dominated by type 1. Mid-west far from center is dominated by type 2 whereas Mid-west closer to center is dominated by type 1. Northwest is dominated by type 1 as well. The center is dominated by type 1 and 2. The Far-Southwest is dominated by type 1.
The next hypthesis is that the feature Sector matches these North-east-south-west areas.
fig, axs = plt.subplots(1, 2,figsize=(20,7))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
for att in ['Sector','CLASE']:
ax = axs[i]
sc = axs[i].scatter(df.Longitude,df.Latitude, s=1, c=df.loc[:,att], cmap=cm)
axs[i].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
axs[i].axis('equal')
fig.colorbar(sc,ax=ax)
i+=1
The sectors corresponds well to the description above, only a slight different division. Sector 1 - Northwest, sector 2 - Mid-west, sector 3 - Northeast etc. The sector indexing matches the houshold grouping well, hence sector can be a good predictor.
The household characteristics of the areas dominated by a household type are now explored. The hypothesis is that high income and high rent areas are related to houshold types in the area.
df_temp=df
fig, axs = plt.subplots(1, 3,figsize=(20,4))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
for att in ['MontoArr','IngresoHogar','CLASE']:
ax = axs[i]
sc = axs[i].scatter(df_temp.Longitude,df_temp.Latitude, s=1, c=df_temp.loc[:,att], cmap=cm)
axs[i].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
axs[i].axis('equal')
fig.colorbar(sc,ax=ax)
i+=1
There is a clear trend in higher rent (MontoArr) and higher income (IngresoHogar) in the Northeast (sector 3). Though, with the much higher values here, than rest of the map it is hard to differentiate. To account for this the lowest rents and incomes are temporarily removed.
#Lowest values are temporarily removed for better coloring
MontoArr_l = df.MontoArr.quantile(0.75)
Ingresso_l = df.IngresoHogar.quantile(0.75)
df_temp = df[(df.MontoArr > MontoArr_l) & (df.MontoArr > MontoArr_l)]
fig, axs = plt.subplots(1, 3,figsize=(20,4))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
for att in ['MontoArr','IngresoHogar','CLASE']:
ax = axs[i]
sc = axs[i].scatter(df_temp.Longitude,df_temp.Latitude, s=1, c=df_temp.loc[:,att], cmap=cm)
axs[i].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
axs[i].axis('equal')
fig.colorbar(sc,ax=ax)
i+=1
When removing the lowest rents and incomes the vast majority of househods are type 3. Hence, it is confirmed that high income and high rent areas are related to houshold types in the area.
The opposite is tested as well.
#Highest values are temporarily removed for better coloring
MontoArr_h = df.MontoArr.median()
Ingresso_h = df.IngresoHogar.median()
df_temp = df[(df.MontoArr < MontoArr_l) & (df.MontoArr < MontoArr_l)]
fig, axs = plt.subplots(1, 3,figsize=(20,4))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
for att in ['MontoArr','IngresoHogar','CLASE']:
ax = axs[i]
sc = axs[i].scatter(df_temp.Longitude,df_temp.Latitude, s=1, c=df_temp.loc[:,att], cmap=cm)
axs[i].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
axs[i].axis('equal')
fig.colorbar(sc,ax=ax)
i+=1
When removing the highest rents and incomes the vast majority of househods are type 1. Hence, it is confirmed that low income and low rent areas are related to houshold types in the area.
IngresoHogar is a household characteristic and therefor not suited for the prediction. Instead Ingreso_Promedio_Zona_MM (Average zonal income (Million Chilean Peso)), can be explaining the income factor as it is clearly correlated with the household type as indicated in the maps below.
fig, axs = plt.subplots(1, 2,figsize=(20,7))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
for att in ['CLASE','Ingreso_Promedio_Zona_MM']:
ax = axs[i]
sc = axs[i].scatter(df.Longitude,df.Latitude, s=1, c=df.loc[:,att], cmap=cm)
axs[i].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
axs[i].axis('equal')
fig.colorbar(sc,ax=ax)
i+=1
The next part of the geographic exploration is to investigate pattern in area charateristics. The first hypothesis is that areas with high levels of crime are central areas.
#Removing NaN replaced coomunes:
crim_mean = df.crimes.mean()
df_c = df[(df.crimes < crim_mean-2) | (df.crimes > crim_mean+2)]
fig, ax = plt.subplots(figsize=(15,7))
cm = plt.cm.get_cmap('RdYlBu_r')
sc = ax.scatter(df_c.Longitude,df_c.Latitude, s=1, c=df_c['crimes'], cmap=cm)
ax.set_title('Crimes by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
fig.colorbar(sc,ax=ax);
Unfortunaly, there is no crime records for the central communes, hence the hypothesis cannot be confirmed nor rejected. An interesting observation is that some outer communens have much higher level of crime than their neighbouring communes. This could e.g. indicate violence in ghettoes or burglary and car thefts in the wealthy neighbourhoods. As this data is not extensive for a throughout geographical investigation, it is conluded here.
Next step is to explore the geographic relationships of the remaining area characteristics.
fig, axs = plt.subplots(4, 2,figsize=(20,20))
cm = plt.cm.get_cmap('RdYlBu_r')
i = 0
j = 0
for att in ['CLASE','Calid_EOD_norm_inv', 'Sup_Prom_Constr_Hab_EOD','DensConstr_EOD','Tiempo_Com_Stgo',
'Dist_est_Metro_MC_KM', 'Dist_salida_Autop_MC_KM', 'Acc_Servicios_tpte_pub']:
ax = axs[i, j]
sc = axs[i, j].scatter(df.Longitude,df.Latitude, s=1, c=df.loc[:,att], cmap=cm)
axs[i, j].set_title('%s by geography' %att)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
fig.colorbar(sc,ax=ax)
if j == 1:
i+=1
j=0
else:
j+=1
Calid_EOD_norm_inv indicates the average construction quality of residential units in the sector, where 0 is highest quality. Areas dominated by household type 3 have highest construction quality. Areas dominated by household type 1 have lowest construction quality.
Sup_Prom_Constr_Hab_EOD inidactes average built surface of residential units in zone. And DensConstr_EOD indicates building density. The wealthy Northeast area is the area with highest density of residential units, but not the most dense area. This indicates that this area is moslty residential, not commercial or industrial. However, the Northeast part closest to the center has a density above 1 indicating towers and several level buildings. The DensConstr_EOD clearly illustrates the center of the city.
Tiempo_Com_Stgo illustrates the time to center. Apart from a few areas (e.g. north of center) the time to the center is dependent on the distance. This is well explained by the distance to the metro stations, Dist_est_Metro_MC_KM, which is evenly distributed accros most of the city. For the area north of the center the connectivity is signifant faster yielding a short time to center even though it is realtive far away. This is an interesting insight for the real estates agency.
Dist_salida_Autop_MC_KM, distance to the highway exit and Acc_Servicios_tpte_pub, accesibility to public transport illustrates a star-like pattern with 5 main routes from the center. This can be compared to the 'the finger-plan' of Copenhagen, see picture, which was presented in 1947 to manage accesibility in a growing city. Areas along the 5 'fingers' have high accesibility to public transports, similar to the plot of accesibility to public transport in Chile. The wealthy area in the Northeast have high accesibility to public transport.
The geographic exploration yields valuable insights for the real estate companies. The wealth charateristics clearly correlates with household types. Hence, Ingreso_Promedio_Zona_MM and MontoArr are promising area characteristic for the prediction.
Calid_EOD_norm_inv, Sup_Prom_Constr_Hab_EOD, DensConstr_EOD are highly promising area characteristics for the prediction as well.
Dist_salida_Autop_MC_KM and Acc_Servicios_tpte_pub might be useful as predictors for the city center, but they are not highly correlated with the type of household according to the geographic exploration.
Dist_est_Metro_MC_KM and Tiempo_Com_Stgo does not seem promising for predictions in city center as these are equal distributed here.
e) Segmentation¶
This section attempts to segment the real estate market based on area characterics. Segmenting the market can help the real estate agency to target specific groups.
The previous four sections have revealed several key insights for the real estate agency, based on data as-is. This section will apply clustering methods to detect data patterns based on multiple features at the same time.
The hypthesis is that the real estate market can be clearly divided into clusters based on area characteristics and that each segment will match some CLASEs more than others. Hence the clusters are created and matched with CLASE.
Segmenting real estate market - standardised data¶
Initially the data is prepared for clustering by One-Hot encoding the Sector variables.
df_e = df
df_e = pd.get_dummies(df, columns = ['Sector'])
print(df_e.shape)
df_e.head(2)
CLASE, Ingresohogar and geographic features are removed and the data is standardised. Initially comunas_id was One-hot encode like the comunas_id, but led to, very long computational time due to a large number of categories. Hence it is removed here too.
from sklearn.preprocessing import StandardScaler
X_e = df_e.drop(['CLASE','IngresoHogar','comunas_id','Longitude','Latitude'],1)
X_e_std = StandardScaler().fit_transform(X_e)
The data is clustered using K-means from sklearn, which clusters data by minimizing the sum of squared distances of observations to K cluster centers.
The approach is inspired by - https://blog.cambridgespark.com/how-to-determine-the-optimal-number-of-clusters-for-k-means-clustering-14f27070048f
from sklearn.cluster import KMeans
Sum_of_squared_distances = []
K = range(1,15)
for k in K:
km_std = KMeans(n_clusters=k,random_state=42)
km_std = km_std.fit(X_e_std)
Sum_of_squared_distances.append(km_std.inertia_)
plt.figure(figsize=(20,10))
plt.plot(K, Sum_of_squared_distances, 'bx-')
plt.xlabel('k')
plt.ylabel('Sum_of_squared_distances')
plt.title('Elbow Method For Optimal k')
plt.show()
By applying the elbow method to select the number of clusters, k=7 clusters is selected, as each of the first 7 clusters reduces the sum of squared distances significantly compared to introducing more clusters.
k=3 clusters would be great for matching with the 3 CLASEs. However for now k=7 seems to be a better option.
The model is fitted on the standardised data:
k=7
km_std=KMeans(k,random_state=42)
km_std.fit(X_e_std)
df['cluster_std']=km_std.fit_predict(X_e_std)
The cluster sizes are illustrated below.
colors = ['b', 'y', 'r','k','c','m','y'] #, 'k', 'w'
plt.figure(figsize=(18,4))
plt.bar(np.arange(k), df.groupby(['cluster_std']).count()['CLASE'].values,color=colors)
plt.xlabel('Cluster')
plt.ylabel('Households')
plt.show()
Below the clusters are linked to their geography.
fig, ax = plt.subplots(figsize=(20,9))
for cluster in range(k):
ax.scatter(df[df.cluster_std==cluster].Longitude,df[df.cluster_std==cluster].Latitude, s=1, c=colors[cluster]);
ax.set_title('Clusters by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal');
The clusters are well divided geographically excatly equal to the sectors. The sectors are removed and the clustering is done all over.
df_e = df
X_e = df_e.drop(['CLASE','IngresoHogar','comunas_id','Longitude','Latitude','Sector'],1)
X_e_std = StandardScaler().fit_transform(X_e)
Sum_of_squared_distances = []
K = range(1,15)
for k in K:
km_std = KMeans(n_clusters=k,random_state=42)
km_std = km_std.fit(X_e_std)
Sum_of_squared_distances.append(km_std.inertia_)
plt.figure(figsize=(20,10))
plt.plot(K, Sum_of_squared_distances, 'bx-')
plt.xlabel('k')
plt.ylabel('Sum_of_squared_distances')
plt.title('Elbow Method For Optimal k')
plt.show()
This time there is no sharp angles, hence the elbow method is not obvious. Applying the elbow method again, k = 7 is selected. k = 2 or 3 would also make sense. After trying each, the k=7 is selected as it seems to have relevance for explaining real estate segmentation patterns.
k=7
km_std=KMeans(k,random_state=42)
km_std.fit(X_e_std)
df['cluster_std']=km_std.fit_predict(X_e_std)
fig, ax = plt.subplots(figsize=(20,9))
for cluster in range(k):
ax.scatter(df[df.cluster_std==cluster].Longitude,df[df.cluster_std==cluster].Latitude, s=1, c=colors[cluster]);
ax.set_title('Clusters by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal');
The clusters are geographiccaly very well divided.
fig, axs = plt.subplots(4, 2,figsize=(20,25))
i = 0
j = 0
for cluster in range(k):
axs[i, j].scatter(df[df.cluster_std==cluster].Longitude,df[df.cluster_std==cluster].Latitude, s=2, c=colors[cluster]);
axs[i, j].set_title('Cluster %i by geography' %cluster)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
if j == 1:
i+=1
j=0
else:
j+=1
Geographiccaly the clusters are very well divided. E.g. one cluster is predominantly located in Northeast Santiago, another cluster is predominantly the Mid- and Northwest suburbs as well as Mid- and SouthEast suburbs. Some clusters are predominantly covering rural districts. Below it is investigated which other characteristics are describing the clusters.
clust_clase_std = np.zeros((3,k))
for i in range(3):
for j in range(k):
clust_clase_std[i,j] = len(df[(df.CLASE==i+1) & (df.cluster_std==j)])
clust_clase_std = clust_clase_std.T
x_std = np.arange(clust_clase_std.shape[0])
dx_std = (np.arange(clust_clase_std.shape[1])-clust_clase_std.shape[1]/3.)/(clust_clase_std.shape[1]+3.)
d_std = 1./(clust_clase_std.shape[1]+3.)
fig, ax = plt.subplots(figsize=(20,9))
for i in range(clust_clase_std.shape[1]):
ax.bar(x_std+dx_std[i],clust_clase_std[:,i], width=d_std, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
plt.title('CLASE per cluster')
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
clust_clase_std_r = clust_clase_std/clust_clase_std.sum(axis=1)[:,None]
fig, ax = plt.subplots(figsize=(20,9))
for i in range(clust_clase_std.shape[1]):
ax.bar(x_std+dx_std[i],clust_clase_std_r[:,i], width=d_std, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
plt.title('CLASE ratio per cluster')
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
The two clusters dominated by CLASE 1 followed up by CLASE 2 are very similar in terms of household types. The three clusters dominated by CLASE 2 followed up by CLASE 1 are very similar in terms of household types too. One cluster is dominated bt CLASE 2 followed up by CLASE 1. Lastly a single clusters is dominated by CLASE 1.
The hypthesis was that the real estate market can be clearly divided into clusters based on area characteristics and that each segment will match one CLASE more than the others is partially confirmed. The clusters are not completely clean, but there are clearly higher ratios of some clases in some clusters than others.
The clusters are now checked for other features to reveal some of these insights. The hypothesis is that each cluster can be described by features of wealth and area characteristics.
Initially wealth features:
colors = ['b', 'y', 'r','k','c','m'] #, 'k', 'w'
for att in ['MontoArr', 'IngresoHogar']:
plt.figure(figsize=(18,4))
plt.bar(np.arange(k), df.groupby(['cluster_std']).mean()[att].values,color=colors)
plt.title('Cluster distribution of - %s' %att)
plt.xlabel('Cluster')
plt.ylabel('Pesos')
plt.show()
The clusters divide income and rent well. E.g. one cluster represents high income and high rent. The relative income and rent is even for all clusters.
The area characterics are checked below.
fig, axs = plt.subplots(5, 2,figsize=(20,20))
i = 0
j = 0
for att in ['Sup_Prom_Constr_Hab_EOD', 'Calid_EOD_norm_inv', 'DensConstr_EOD',
'Dist_est_Metro_MC_KM', 'Dist_salida_Autop_MC_KM', 'Tiempo_Com_Stgo','Ingreso_Promedio_Zona_MM',
'Acc_Servicios_tpte_pub','crimes']:
axs[i,j].bar(np.arange(k), df.groupby(['cluster_std']).mean()[att].values,color=colors)
axs[i,j].set_title('Cluster average - %s' %att)
if j == 1:
i+=1
j=0
else:
j+=1
The clustering reveals patterns very similar to what the patterns for the majority CLASE. The clusters can be described according to their features values, and the real estate agency can use this to find real estate that matches a household knowing its income, rent, prefered distance to metro etc.
E.g. one cluster is mixed of all clases and has mid-low rent, mid-low income, half industrial/commercial area half resdiential area, low construction quality, high density, medium distance to metro, low distance to highway and low average income in the neighbourhood.
The hypothesis was that each cluster can be described by features of wealth and area characteristics and is confirmed. If a real estate agent gets a real estate with the area characteristics matching one of the clusters, he can send an advertisement to households matching this cluster. This way of targeting will improve the efforts of the real estate agents.
Segmenting real estate market - PCA transformed data¶
Here it is checked whether using principal components for clustering will yield a more clear distinction of househould type by clusters than the previous section. The data is initially PCA transformed.
from sklearn import decomposition
pca=decomposition.PCA()
pca.fit(X_e)
expl=pca.explained_variance_ratio_
cdf=[sum(expl[:i+1]) for i in range(len(expl))]
plt.figure(figsize=(18,4))
plt.plot(range(len(expl)), cdf);
for i in range(4):
print('PC-%i explains: %.5f of the variance' % (i, cdf[i]))
The first principal component explains 96.5% of the variance and together with the second PC they explain 99.7% of the variance. The PCs are used for clustering with only 2 dimensions, the clustering can now be illustated.
X_pca=pca.transform(X_e)
plt.figure(figsize=(18,6))
plt.scatter(X_pca[:, 0],X_pca[:,1]);
The scatter doesn't reveal any obvioues clusters apart from the far right low density area. The transformed data is clustered using KMeans as before.
from sklearn.cluster import KMeans
Sum_of_squared_distances = []
K = range(1,15)
for k in K:
km_pca = KMeans(n_clusters=k,random_state=42)
km_pac = km_pca.fit(X_pca[:,0:2])
Sum_of_squared_distances.append(km_pca.inertia_)
plt.figure(figsize=(20,10))
plt.plot(K, Sum_of_squared_distances, 'bx-')
plt.xlabel('k')
plt.ylabel('Sum_of_squared_distances')
plt.title('Elbow Method For Optimal k')
plt.show()
By applying the elbow method to select the number of clusters, k=3 clusters is selected., though k=2 or 4 are also promising options.
k=3
km_pca=KMeans(k,random_state=42)
km_pca.fit(X_pca)
df['cluster_pca']=km_pca.fit_predict(X_pca)
df['X_PC1'] = X_pca[:,0]
df['X_PC2'] = X_pca[:,1]
colors = ['m', 'g', 'y'] #,'k','c','m','y'] #, 'k', 'w'
plt.figure(figsize=(18,4))
plt.bar(np.arange(k), df.groupby(['cluster_pca']).count()['CLASE'].values,color=colors)
plt.xlabel('Cluster')
plt.ylabel('Households')
plt.show()
plt.figure(figsize=(18,6))
plt.scatter(df.X_PC1,df.X_PC2,c=df.cluster_pca);
The 3 clusters are highly imbalanced. The scatter indicates some fair patterns to differentitate.
fig, ax = plt.subplots(figsize=(20,9))
for cluster in range(k):
ax.scatter(df[df.cluster_pca==cluster].Longitude,df[df.cluster_pca==cluster].Latitude, s=1, c=colors[cluster]);
ax.set_title('Clusters by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal');
The geogography is again clearly divided by the clusters.
clust_clase_pca = np.zeros((3,k))
for i in range(3):
for j in range(k):
clust_clase_pca[i,j] = len(df[(df.CLASE==i+1) & (df.cluster_pca==j)])
clust_clase_pca = clust_clase_pca.T
clust_clase_pca_r = clust_clase_pca/clust_clase_pca.sum(axis=1)[:,None]
x_pca = np.arange(clust_clase_pca.shape[0])
dx_pca = (np.arange(clust_clase_pca.shape[1])-clust_clase_pca.shape[1]/3.)/(clust_clase_pca.shape[1]+3.)
d_pca = 1./(clust_clase_pca.shape[1]+3.)
fig, ax = plt.subplots(figsize=(20,9))
for i in range(clust_clase_pca.shape[1]):
ax.bar(x_pca+dx_pca[i],clust_clase_pca_r[:,i], width=d_pca, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
plt.title('CLASE ratio per cluster')
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
This clustering can been interpreted as one cluster being the upper class, one clsuter being the upper middle class and one cluster the lower middle class.
The two clustering approaces reveal valuable insights for the real estate agency. Neither is perfectly dividing the CLASE categories, but each clustering appraoch segments the market in categories driven by different level of area characteristics and household charateristics.
The hypothesis was that this PCA approach would yield better distinction in cluster with regards to CLASE, this is confirmed.
Segmenting real estate market - standardised and PCA transformed data¶
Lastly the data is standardised and PCA transformed before clustering, to see what impact this will have on the business value for the real estate company. The hypothesis is that clustering on the standardised and PCA transformed data will lead to higher purity of clusters with regards to CLASE.
pcastd=decomposition.PCA()
pcastd.fit(X_e_std)
expl=pcastd.explained_variance_ratio_
cdf=[sum(expl[:i+1]) for i in range(len(expl))]
plt.figure(figsize=(18,4))
plt.plot(range(len(expl)), cdf);
for i in range(11):
print('PC-%i explains: %.5f of the variance' % (i, cdf[i]))
PC-0 to PC-6 is explaining 90% of the variance, hence these are applied. To illustrate the three first PCs are used for a 3D plot below.
from mpl_toolkits.mplot3d import Axes3D
X_pcastd=pcastd.transform(X_e_std)
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_pcastd[:, 0],X_pcastd[:,1],X_pcastd[:,2]);
ax.set_xlabel('PC0')
ax.set_ylabel('PC1')
ax.set_zlabel('PC2')
ax.view_init(15, 90)
plt.show()
Sum_of_squared_distances_pcastd = []
K = range(1,15)
for k in K:
km_pcastd = KMeans(n_clusters=k)
km_pcastd = km_pcastd.fit(X_pcastd[:,0:7]) #PC-0 to PC-6
Sum_of_squared_distances_pcastd.append(km_pcastd.inertia_)
plt.figure(figsize=(20,10))
plt.plot(K, Sum_of_squared_distances_pcastd, 'bx-')
plt.xlabel('k')
plt.ylabel('Sum_of_squared_distances')
plt.title('Elbow Method For Optimal k')
plt.show()
There is no clear elbow to select the optimal k, but after k=3 the curve flattens somewhat hence k=3 is selected.
k=3
km_pcastd=KMeans(k,random_state=42)
km_pcastd.fit(X_pcastd)
df['cluster_pcastd']=km_pcastd.fit_predict(X_pcastd)
plt.figure(figsize=(18,4))
plt.bar(np.arange(k), df.groupby(['cluster_pcastd']).count()['CLASE'].values,color=colors)
plt.xlabel('Cluster')
plt.ylabel('Households')
plt.show()
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_pcastd[:, 0],X_pcastd[:,1],X_pcastd[:,2],c=df.cluster_pcastd);
ax.set_xlabel('PC0')
ax.set_ylabel('PC1')
ax.set_zlabel('PC2')
ax.view_init(15, 90)
plt.show()
The clusters seem to be well seperated by the first three PCs. The geographic link is illustrated below:
fig, ax = plt.subplots(figsize=(20,9))
for cluster in range(k):
ax.scatter(df[df.cluster_pcastd==cluster].Longitude,df[df.cluster_pcastd==cluster].Latitude, s=1, c=colors[cluster]);
ax.set_title('Clusters by geography')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal');
Geographically the clusters are very will divided ass well. Lastly the clusters are mathced to CLASE.
clust_clase_pcastd = np.zeros((3,k))
for i in range(3):
for j in range(k):
clust_clase_pcastd[i,j] = len(df[(df.CLASE==i+1) & (df.cluster_pcastd==j)])
clust_clase_pcastd = clust_clase_pcastd.T
clust_clase_pcastd_r = clust_clase_pcastd/clust_clase_pcastd.sum(axis=1)[:,None]
x_pcastd = np.arange(clust_clase_pcastd.shape[0])
dx_pcastd = (np.arange(clust_clase_pcastd.shape[1])-clust_clase_pcastd.shape[1]/3.)/(clust_clase_pcastd.shape[1]+3.)
d_pcastd = 1./(clust_clase_pcastd.shape[1]+3.)
fig, ax = plt.subplots(figsize=(20,9))
for i in range(clust_clase_pcastd.shape[1]):
ax.bar(x_pcastd+dx_pcastd[i],clust_clase_pcastd_r[:,i], width=d_pcastd, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
plt.title('CLASE ratio per cluster')
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
The combination of standardisation and PCA tranformation yields clusters that divides household types well too. In order to confirm the hypthosesis of this clustering being better, the last part of this analysis will compare the three methods on purity with regards to CLASE.
fig, ax = plt.subplots(3, 1,figsize=(20,20))
for i in range(clust_clase_std.shape[1]):
ax[0].bar(x_std+dx_std[i],clust_clase_std_r[:,i], width=d_std, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[0].set_title('std - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_std_r,1))/sum(sum(clust_clase_std_r))))
for i in range(clust_clase_pca.shape[1]):
ax[1].bar(x_pca+dx_pca[i],clust_clase_pca_r[:,i], width=d_pca, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[1].set_title('pca - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_pca_r,1))/sum(sum(clust_clase_pca_r))))
for i in range(clust_clase_pcastd.shape[1]):
ax[2].bar(x_pcastd+dx_pcastd[i],clust_clase_pcastd_r[:,i], width=d_pcastd, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[2].set_title('stdpca - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_pcastd_r,1))/sum(sum(clust_clase_pcastd_r))))
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
According to this result, PCA without standardisation yielded highest purity. For the standardisation only model, it was also deemed reasonable to apply k=3. This will be done to have a more fair comparison.
k=3
km_std3=KMeans(k,random_state=42)
km_std3.fit(X_e_std)
df['cluster_std3']=km_std3.fit_predict(X_e_std)
#matrix
clust_clase_std3 = np.zeros((3,k))
for i in range(3):
for j in range(k):
clust_clase_std3[i,j] = len(df[(df.CLASE==i+1) & (df.cluster_std3==j)])
clust_clase_std3 = clust_clase_std3.T
clust_clase_std3_r = clust_clase_std3/clust_clase_std3.sum(axis=1)[:,None]
fig, ax = plt.subplots(3, 1,figsize=(20,20))
for i in range(clust_clase_std3.shape[1]): #Obs x_pca, dx_pca and d_pca have the right dims here
ax[0].bar(x_pca+dx_pca[i],clust_clase_std3_r[:,i], width=d_pca, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[0].set_title('std - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_std3_r,1))/sum(sum(clust_clase_std3_r))))
for i in range(clust_clase_pca.shape[1]):
ax[1].bar(x_pca+dx_pca[i],clust_clase_pca_r[:,i], width=d_pca, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[1].set_title('pca - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_pca_r,1))/sum(sum(clust_clase_pca_r))))
for i in range(clust_clase_pcastd.shape[1]):
ax[2].bar(x_pcastd+dx_pcastd[i],clust_clase_pcastd_r[:,i], width=d_pcastd, color= ['xkcd:sky blue','xkcd:dark blue','xkcd:bright blue'][i], label="{}".format(['CLASE 1','CLASE 2','CLASE 3'][i]))
ax[2].set_title('stdpca - CLASE ratio per cluster. Purity score = %.2f' %(sum(np.max(clust_clase_pcastd_r,1))/sum(sum(clust_clase_pcastd_r))))
plt.xlabel('Cluster')
plt.legend(framealpha=1)
plt.show()
With this change the Kmean(3) on standardised data and the Kmean(3) on standardised and PCA transformed data clusters equally. The Kmean(3) on non standardised data but PCA transformed yield the best purity by being significantly better at clustering CLASE 3 by itself. This can be a result of the very high value for neighbourhood income and rent having higher influence whitout standardisation, and these features are primarliy the drivers for CLASE 3. However the standardised models yields clusters that are dominated by each CLASE which might be more valuable for the real estate agency.
The three key take-aways from this analysis are:
1) Different segmentations the market are established for the real estate agency to target the right buyers.
2) The PCA clustering without standardisation is better at segmenting CLASE 3.
3) The standardised clusterings are better at segmenting with a segment dominated by each house hold type.
Since this clustering is based on the entire data of the city it will not be used to predict in the next section. The model could be trained for the prediction challenge to, but other machine learning methods are applied there to widen the toolbox.
Prediction challenge¶
The objective is to predict the type of household that will buy a real estate on basis of its area characteristics. This prediction can help real estate companies to target their advertising to the segment with highest potential.
To evaluate models, they are build as multicategorical calssifications, predicting the most likely category of household type to buy a real estate. However, the real estate agencies might be interested in the second and third most likely household types as well, if the probability of these are high. Therefore the output of the final model will be in the format of:
House X is relevant for households:
Type 1 - y%
Type 2 - z%
Type 3 - w%
The model will be trained on data outside the center of Santiago, and tested on the center of Santiago. The center is defined as the area within the coordinates 70.7E, 33.4S to 70.6E, 33.5S. The categorical feature Sector is one-hot encoded. To illustrate the train and test split the geometry is merged back to the data set. The dataframe used here is the df_all dataframe as it has all features and has been cleaned.
df_all = df_all.reset_index()
df_all = df_all.merge(coords_and_geom[['Hogar','geometry']], how='left',on='Hogar')
df_all = df_all.set_index('Hogar')
df_all = pd.get_dummies(df_all, columns = ['Sector'])
df_all.head()
For predicting the target feature, 'CLASE', three different set of features will be made.
1) Test and train data where the 'acc' features found in the correlation matrix are removed.
2) Test and train data where the 'acc' features are represented by PCA vectores.
3) Test and train data where all features are presented by PCA vectors.
Before the data is split the PCA vectores are made and added to the dataframe. First the PCA vectors for the 'acc' features are made.
# Acc feature names
acc_feature_names = ['Acc_Comercio_tpte_pub','Acc_Educacion_tpte_pub', 'Acc_Habitacion_tpte_pub',
'Acc_Industria_tpte_pub', 'Acc_Servicios_tpte_pub', 'Acc_Comercio_auto',
'Acc_Educacion_auto', 'Acc_Habitacion_auto', 'Acc_Industria_auto','Acc_Servicios_auto',]
# New data frame only with the acc features
acc_features = df_all[acc_feature_names]
# Standardize the data
acc_features_std = StandardScaler().fit_transform(acc_features)
# Perform PCA:
pca = decomposition.PCA()
pca.fit(acc_features_std)
expl=pca.explained_variance_ratio_
cdf=[sum(expl[:i+1]) for i in range(len(expl))]
plt.figure(figsize=(18,4))
plt.plot(range(len(expl)), cdf);
for i in range(6):
print('PC-%i explains: %.5f of the variance' % (i+1, cdf[i]))
The first 4 principal components are capable of explaining 99.6% of the variance of the 'acc' features. These 4 vectors are therefor used to explain the 'acc' features.
acc_pca = pca.transform(acc_features_std)
df_all['PCA1_acc'] = acc_pca[:,0]
df_all['PCA2_acc'] = acc_pca[:,1]
df_all['PCA3_acc'] = acc_pca[:,3]
df_all['PCA4_acc'] = acc_pca[:,4]
A PCA for all the features in the set is now made.
# Feature names to use for PCA
features = ['Sup_Prom_Constr_Hab_EOD','MontoArr',
'Calid_EOD_norm_inv', 'DensConstr_EOD', 'Dist_est_Metro_MC_KM',
'Dist_salida_Autop_MC_KM', 'Tiempo_Com_Stgo',
'Ingreso_Promedio_Zona_MM', 'Acc_Comercio_tpte_pub',
'Acc_Educacion_tpte_pub', 'Acc_Habitacion_tpte_pub',
'Acc_Industria_tpte_pub', 'Acc_Servicios_tpte_pub', 'Acc_Comercio_auto',
'Acc_Educacion_auto', 'Acc_Habitacion_auto', 'Acc_Industria_auto',
'Acc_Servicios_auto', 'Latitude', 'Longitude', 'crimes',
'Sector_1', 'Sector_2', 'Sector_3', 'Sector_4', 'Sector_5',
'Sector_6', 'Sector_7']
# New data frame only with the features choosen above
df_features = df_all[features]
# Standardize the data
df_features_std = StandardScaler().fit_transform(df_features)
# Perform PCA:
pca = decomposition.PCA()
pca.fit(df_features_std)
expl=pca.explained_variance_ratio_
cdf=[sum(expl[:i+1]) for i in range(len(expl))]
plt.figure(figsize=(18,4))
plt.plot(range(len(expl)), cdf);
for i in range(8,16):
print('PC-%i explains: %.5f of the variance' % (i, cdf[i]))
We see that more than 95% of the variance is explained by the first 11 PCA components. These features are also saved to the dataframe.
all_pca = pca.transform(df_features_std)
# Define alle PCA's
for i in range(0,11):
df_all['PCA-%i' % (i+1)] = all_pca[:,i]
df_all.head()
With all PCA vectors in the dataframe the data can now be divided into train and test data. It should be noted that 'IngresoHogar'is removed since it is a household characteristic which is not a real estate characteristic.
#Defining independent variables X, and dependent variables y. (Longitude and Latitude are temporarily kept to map below)
X_att = ['Sup_Prom_Constr_Hab_EOD','MontoArr',
'Calid_EOD_norm_inv', 'DensConstr_EOD', 'Dist_est_Metro_MC_KM',
'Dist_salida_Autop_MC_KM', 'Tiempo_Com_Stgo',
'Ingreso_Promedio_Zona_MM', 'Acc_Comercio_tpte_pub',
'Acc_Educacion_tpte_pub', 'Acc_Habitacion_tpte_pub',
'Acc_Industria_tpte_pub', 'Acc_Servicios_tpte_pub', 'Acc_Comercio_auto',
'Acc_Educacion_auto', 'Acc_Habitacion_auto', 'Acc_Industria_auto',
'Acc_Servicios_auto', 'Latitude', 'Longitude', 'crimes','geometry',
'Sector_1', 'Sector_2', 'Sector_3', 'Sector_4', 'Sector_5',
'Sector_6', 'Sector_7', 'PCA1_acc','PCA2_acc','PCA3_acc','PCA4_acc',
'PCA-1','PCA-2','PCA-3','PCA-4','PCA-5','PCA-6','PCA-7','PCA-8','PCA-9','PCA-10','PCA-11']
y_att = ['CLASE','Longitude','Latitude']
# ------- Train and Test split ----------- #
#Define borders
TestLongMin = -70.7
TestLongMax = -70.6
TestLatMin = -33.5
TestLatMax = -33.4
# Split sets
X_train = df_all[(df_all.Longitude < TestLongMin) | (TestLongMax < df_all.Longitude) | (df_all.Latitude < TestLatMin) | (TestLatMax < df_all.Latitude)][X_att]
X_test = df_all[(df_all.Longitude >= TestLongMin) & (TestLongMax >= df_all.Longitude) & (df_all.Latitude >= TestLatMin) & (TestLatMax >= df_all.Latitude)][X_att]
y_train = df_all[(df_all.Longitude < TestLongMin) | (TestLongMax < df_all.Longitude) | (df_all.Latitude < TestLatMin) | (TestLatMax < df_all.Latitude)][y_att]
y_test = df_all[(df_all.Longitude >= TestLongMin) & (TestLongMax >= df_all.Longitude) & (df_all.Latitude >= TestLatMin) & (TestLatMax >= df_all.Latitude)][y_att]
#Shapes of train and test sets
print('X_train shape:',X_train.shape)
print('X_test shape:',X_test.shape)
print('y_train shape:',y_train.shape)
print('y_test shape:',y_test.shape)
print('%.1f of the data is train data.' %(X_train.shape[0]/(X_train.shape[0]+X_test.shape[0])*100))
#Define geodataframe and plot
df_geo = gdp.GeoDataFrame(df_all, crs = utm, geometry = df_all.geometry)
plt.figure(figsize=(15,10));
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
plt.scatter(X_train.Longitude,y_train.Latitude, s=1, c='tab:blue')
plt.scatter(X_test.Longitude,y_test.Latitude, s=1, c='red')
plt.title('Coordinates of real estate transaction data in Santiage de Chile - Training data (blue) and test data (red)');
plt.show()
#<W ^N E>
The latitude, longitude and geometry columns are removed from the y as we only want CLASE. Furthermore the longitude and latitude for the test set is saved for later use.
test_coor = y_test.loc[:,'Longitude':'Latitude']
y_train = y_train.drop(['Longitude','Latitude'],1)
y_test = y_test.drop(['Longitude','Latitude'],1)
As stated in the initial data analysis, the target feature was imbalanced. Here the distribution is checked for both the training and test set.
p1 = y_train.reset_index()
p2 = y_test.reset_index()
## Set subplot
fig, axes = plt.subplots(1, 2, figsize=(15,5))
axes[0].set_title('Training targets')
p1.groupby('CLASE').count().plot(kind='bar',ax=axes[0])
axes[1].set_title('Test targets')
p2.groupby('CLASE').count().plot(kind='bar',ax=axes[1])
plt.show()
The training data is highly based on CLASE 1 and 2, whereas the test data balanced towards CLASE 2 and 3. The train and test set are fixed so the distribution cannot be changed. This can potentially cause problems as our model will be trained on one distribution of the data but it will be evaluated on a different distribution. During the prediction the scoring metric for the models will be changed to take this into account.
The 3 described set of features are now defined.
1) Test and train data where the 'acc' features found in the correlation matrix are removed.
X1_features = ['Sup_Prom_Constr_Hab_EOD', 'Calid_EOD_norm_inv', 'DensConstr_EOD',
'Dist_est_Metro_MC_KM', 'Dist_salida_Autop_MC_KM', 'Tiempo_Com_Stgo',
'Ingreso_Promedio_Zona_MM', 'Acc_Servicios_tpte_pub',
'Acc_Servicios_auto', 'crimes', 'Sector_1', 'Sector_2', 'Sector_3',
'Sector_5', 'Sector_6', 'Sector_7','MontoArr']
X1_train = X_train[X1_features]
X1_test = X_test[X1_features]
print(X1_train.shape)
print(X1_test.shape)
The X1 test and train set thus has 17 features.
2) Test and train data where the 'acc' features are represented by PCA vectores.
X2_features = ['Sup_Prom_Constr_Hab_EOD', 'Calid_EOD_norm_inv', 'DensConstr_EOD',
'Dist_est_Metro_MC_KM', 'Dist_salida_Autop_MC_KM', 'Tiempo_Com_Stgo',
'Ingreso_Promedio_Zona_MM','crimes', 'Sector_1', 'Sector_2', 'Sector_3',
'Sector_5', 'Sector_6', 'Sector_7', 'PCA1_acc', 'PCA2_acc', 'PCA3_acc',
'PCA4_acc','MontoArr']
X2_train = X_train[X2_features]
X2_test = X_test[X2_features]
print(X2_train.shape)
print(X2_test.shape)
The X2 test and train set thus has 19 features.
3) Test and train data where all features are presented by PCA vectors.
X3_features = ['PCA-1', 'PCA-2', 'PCA-3', 'PCA-4', 'PCA-5',
'PCA-6', 'PCA-7', 'PCA-8', 'PCA-9', 'PCA-10','PCA-11']
X3_train = X_train[X3_features]
X3_test = X_test[X3_features]
print(X3_train.shape)
print(X3_test.shape)
This leaves 11 features in the X3 train and test set.
The data is standardized before used for modelling. This is done to ensure that we are able to compare features with different units. For X1 all the data is standardized, for X2 all features but the PCA features are standardized and for X3 no features are standardized as this has all ready been done.
# Standardize X1
X1_mean = X1_train.mean(axis=0)
X1_std = X1_train.std(axis=0)
X1_train = (X1_train - X1_mean) / X1_std
X1_test = (X1_test - X1_mean) / X1_std
# Standardize X2
X2_mean = X2_train.iloc[:,0:14].mean(axis=0)
X2_std = X2_train.iloc[:,0:14].std(axis=0)
X2_train.iloc[:,0:14] = (X1_train - X1_mean) / X1_std
X2_test.iloc[:,0:14] = (X1_test - X1_mean) / X1_std
The training data is ready for modelling of classifiers and the test data is ready for evaluation of the classifiers. The two evaluation functions below are used.
The original function originates from the lectures, but is tailored to this project and split in two defintion. confusion_mat plot a confusion matrix while the definition evaluate_scores gives the accuracy, precision, recall and f1 score.
from sklearn.metrics import confusion_matrix, recall_score, accuracy_score, precision_score,classification_report
from exercise_utils import plot_confusion_matrix
# function to evaluate predictions
def confusion_ma(y_true, y_pred):
# calculate and display confusion matrix
labels = np.unique(y_true)
cm = confusion_matrix(y_true, y_pred, labels=labels)
plot_confusion_matrix(y_true, y_pred,classes=labels,normalize=True)
return plt.show()
# function to evaluate only accuracy, precision, recall and f1 score
def evaluate_scores(y_true, y_pred):
# calculate confusion matrix
labels = np.unique(y_true)
cm = confusion_matrix(y_true, y_pred, labels=labels)
# calculate precision, recall, and F1 score
accuracy = float(np.trace(cm)) / np.sum(cm)
precision = precision_score(y_true, y_pred, average=None, labels=labels)[1]
recall = recall_score(y_true, y_pred, average=None, labels=labels)[1]
f1 = 2 * precision * recall / (precision + recall)
print("accuracy:", accuracy)
print("precision:", precision)
print("recall:", recall)
print("f1 score:", f1)
return accuracy, precision, recall, f1
3 classificaiton models will be evaluated:
- Logistics Regression
- K-Nearest Neighbours
- Random Forest
10-fold cross validation (CV) is applied to select the best model. 10-fold CV yields 14941 observations for training and 2929 observations for validation. 10-fold CV is selected as a trade off between the precise but computational heavy Leave-One-Out method and the less precise but computational light Hold-Out method.
Each model will undergo parameter tuning to achieve the best possible version of the model. The parameter ranges are selected by the authors, and will be changed if the best setting is at the edge of the range. The parameter tunning is perfomed using grid search, entailing that all parameter combinations are tested. This appraoch is also called exhaustive parameter tuning.
f1_micro is selected as the scoring metric for the parameter tuning. It gives each sample-class pair an equal contribution to the overall metric. This way it is considered equally important to classify each type of houshold, even though some types are dominating. According to sckikit-learn: Micro-averaging may be preferred in multilabel settings, including multiclass classification where a majority class is to be ignored.. No classes will be ignored with the scoring, but the performance for the largest category CLASE 2, won't count more than the performance of the other categories.
GridSearchCV from sklearn combine the crossvalidation and parameter tunning.
#GridSearch
from sklearn.model_selection import GridSearchCV
#10-fold CV
kf = 10
#Scoring - https://scikit-learn.org/stable/modules/model_evaluation.html
scor = 'f1_micro'
Inspiration for the GridSearchCV and the scoring is obtained here:
https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html
https://scikit-learn.org/stable/modules/model_evaluation.html#scoring-parameter
Scoring types
https://datascience.stackexchange.com/questions/15989/micro-average-vs-macro-average-performance-in-a-multiclass-classification-settin
The sklearn dummy classifier is used as baseline for the model selection process. In the classifier the strategy ='most_frequent' is set as the model then will predict everything as the most frequent label in the training set. It therefor does not matter which X is given as the results only depend on the the y.
from sklearn.dummy import DummyClassifier
#Estimator
dummy_clf = DummyClassifier(strategy = 'most_frequent')
dummy_clf.fit(X1_train,y_train)
#Results
y_pred = dummy_clf.predict(X1_test)
evaluate_scores(y_test,y_pred)
print(classification_report(y_test, y_pred))
confusion_ma(y_test,y_pred)
The model, as it was told, classifies everything as being CLASE = 2 which gives a f1 score of 0.60 which we later will see is quit high. However, if we think from the companys point of view, the f1 score is not nessasary a very good way to score the models. The precision and accuracy for the different classes might be more important. This is especially important for the customers in the CLASE = 3, as these customers, as shown, are customers with high wages and thereby higher income for the company. We therefor want the accuracy to be as high for this class as possible.
With this in mind the DummyClassifier does very badly as it does not classify any of the CLASE = 3 correctly. Neither does it predict any of the CLASE = 1 correctly. If we went ahead with this model we would there for send a lot of wrong porposals out to the CLASE's.
The parameters which are choosen for tuning in the LogisticRegression are: penalty, C and fit_intercept. The penalty specifies which regularization to use. The model is given the l1 and l2 to choose between. The C values specifies the inverse of the regularization strength, lower values of C constraint the model more. Finally the fit_intercept specifies whether or not a constant should be added to the function.
In the logistic regression the following things are further specified:
multi_class = 'ovr'- is specified to avoid error messages.'ovr'is the default setting.solver = 'liblinear'- is specified to avoid error messages.'liblinear'is the default setting.random_state = 42- is specified to make sure the same result is achieved every time.class_weight = 'balanced'- is specified to acommodate the imbalance in the training set.
A logistic regression is now fitted for X1_train.
from sklearn.linear_model import LogisticRegression
#Estimator
logreg = LogisticRegression(multi_class = 'ovr',solver = 'liblinear', class_weight = 'balanced', random_state = 42)
#Parameters for the tuning
logreg_pars = {'penalty': ['l1','l2'],'C' : [0.0009,0.001,0.002], 'fit_intercept': ['True','False']}
#GridSearchCV
logreg_clf1 = GridSearchCV(logreg, logreg_pars, cv=kf, scoring=scor)
logreg_clf1.fit(X1_train,np.ravel(y_train))
#Results
print('Best parameters:',logreg_clf1.best_params_)
print('Best score:',logreg_clf1.best_score_)
A logistic regression is now fitted for X2_train.
#Parameters for the tuning
logreg_pars = {'penalty': ['l1','l2'],'C' : [0.0009,0.001,0.0015], 'fit_intercept': ['True','False']}
#GridSearchCV
logreg_clf2 = GridSearchCV(logreg, logreg_pars, cv=kf, scoring=scor)
logreg_clf2.fit(X2_train,np.ravel(y_train))
#Results
print('Best parameters:',logreg_clf2.best_params_)
print('Best score:',logreg_clf2.best_score_)
Finally a logistic regression is fitted for X3_train.
#Parameters for the tuning
logreg_pars = {'penalty': ['l1','l2'],'C' : [0.027,0.028,0.029], 'fit_intercept': ['True','False']}
#GridSearchCV
logreg_clf3 = GridSearchCV(logreg, logreg_pars, cv=kf, scoring=scor)
logreg_clf3.fit(X3_train,np.ravel(y_train))
#Results
print('Best parameters:',logreg_clf3.best_params_)
print('Best score:',logreg_clf3.best_score_)
From the results of the models all three models choose the to have a constant added to the decision function and that the penalty should be l1. The regularization strength for the three models is different as the two first models choose a much lower C value and thereby a stronger regularization.
With the trained three logistic regressions it is now time to evaluate their performance. We will start of by looking at the models scores.
y_pred1 = logreg_clf1.predict(X1_test)
print('---------- X1 -----------')
logreg_eval1 = evaluate_scores(y_test,y_pred1)
print(classification_report(y_test, y_pred1))
y_pred2 = logreg_clf2.predict(X2_test)
print('\n---------- X2 -----------')
logreg_eval2 = evaluate_scores(y_test,y_pred2)
print(classification_report(y_test, y_pred2))
y_pred3 = logreg_clf3.predict(X3_test)
print('\n---------- X3 -----------')
logreg_eval3 = evaluate_scores(y_test,y_pred3)
print(classification_report(y_test, y_pred3))
From the scores presented above we see that all models have worse f1 score than the dummy classifier. However, all models have a higher overall accuracy and higher overall precision and for each CLASE. Therefor any of these logisticregression models would be much more desirable than the DummyClassifier. We futher see that the model trained on X1 has the highest precision for CLASE = 3 but also that it has the lowest overall accuracy and precision and a low accuracy for CLASE = 3. To be able to decide which model and thereby which dataset is gives the best results with the logisticregression we also need to inspect the confusion matrices.
# Confusion matrix
print('---------- X1 -----------')
confusion_ma(y_test,y_pred1)
print('---------- X2 -----------')
confusion_ma(y_test,y_pred2)
print('---------- X3 -----------')
confusion_ma(y_test,y_pred3)
By inspecting the confusion matrix for X1 we see that eventhough CLASE = 3 had the highest precision this is because it generally just predicts much fewer observations as CLASE = 3 and also has the lowest amount of true positives for this CLASE. The data for X1 is therefor not choosen for the logistic regression. By in comparing the f1 score for CLASE = 3 for the different models the model trained on X1 also has the lowest score.
With regards to the models trained on X2 and X3 the overall accuracy and precision were very similar. By inspecting the two confusion matrices this is verified. We further see that the model trained on X2 only predicts 68% of CLASE = 3 correctly while the model trained on X3 predicts 79% correctly. However the model trained on X2 is very good at not predicting many households which are of CLASE = 1 as CLASE = 3 (only 10%) and the same with households which are CLASE = 3 as CLASE = 1 (4%). This is very important as the people who will buy a household of the type CLASE = 3 might also buy a house of type CLASE = 2, but they are not likely to want a house of type CLASE = 1. For the people who will buy a household of type CLASE = 1 it is the other way around. As for the people who will will buy a household of CLASE = 2 it is more likely that they might buy a household of any of the categories and it is there for okay that larger shares for this CLASE are predicted falsly.
The model which is trained on the X3 data (the data where all vectors are transformed into PCA vectors) is choosen as the best data for the logistic regression. The model has a l1 penalty of 0.028 and uses a constant in the decision function. It predicts the largest share of CLASE = 3 corretly while still predicting a large share of the CLASE = 1 truly. It further has a low level CLASE = 3 which are predicted as CLASE = 1 and CLASE = 1 which are predicted as CLASE = 3. The models also has the highest f1 score for CLASE = 3.
K-Nearest Neighbours
In this section we will look at the K-Nearest Neighbour classifier. For parameter tuning the amount of neighbors and the weights are specified. For the weights the model can choose between uniform which is where points are weighted equally and distance where closer neighbors will have a bigger influence than neighbors which are further away.
For deciding on the amount of neighbors the models are initially trained with 1, 10, 100 and 1000 neighbors to get an idea of what amount of neighbors are needed. With this result a more precise parameter tuning is done after.
A KNN is fitted for X1_train.
from sklearn.neighbors import KNeighborsClassifier
#Estimator
KNclas = KNeighborsClassifier()
#Parameters for the tuning
#KNclas_pars = {'n_neighbors' : [1,10,100,1000], 'weights': ['uniform','distance']}
KNclas_pars = {'n_neighbors' : [1050,1100,1150], 'weights': ['uniform','distance']}
#GridSearchCV
KNclas_clf1 = GridSearchCV(KNclas, KNclas_pars, cv=kf, scoring=scor)
KNclas_clf1.fit(X1_train,np.ravel(y_train))
#Results
print('Best parameters:',KNclas_clf1.best_params_)
print('Best score:',KNclas_clf1.best_score_)
A KNN is now fitted for X2_train
#Parameters for the tuning
#KNclas_pars = {'n_neighbors' : [1,10,100,1000], 'weights': ['uniform']}
KNclas_pars = {'n_neighbors' : [1450,1500,1550], 'weights': ['uniform']}
#GridSearchCV
KNclas_clf2= GridSearchCV(KNclas, KNclas_pars, cv=kf, scoring=scor)
KNclas_clf2.fit(X2_train,np.ravel(y_train))
#Results
print('Best parameters:',KNclas_clf2.best_params_)
print('Best score:',KNclas_clf2.best_score_)
Finally a KNN is fitted for X3_train.
#Parameters for the tuning
#KNclas_pars = {'n_neighbors' : [1,10,100,1000], 'weights': ['uniform','distance']}
KNclas_pars = {'n_neighbors' : [1150,1200,1250], 'weights': ['uniform']}
#GridSearchCV
KNclas_clf3 = GridSearchCV(KNclas, KNclas_pars, cv=kf, scoring=scor)
KNclas_clf3.fit(X3_train,np.ravel(y_train))
#Results
print('Best parameters:',KNclas_clf3.best_params_)
print('Best score:',KNclas_clf3.best_score_)
The three models choose different amount of neighbors as the optimal, however they all choose above 1000 neighbors, all three models also thoose the weights to be uniform.
With the three different KNN fitted it is time to evaluate their performance. We will start by inspecting the model scores.
y_pred1 = KNclas_clf1.predict(X1_test)
print('---------- X1 -----------')
knn_eval1 = evaluate_scores(y_test,y_pred1)
print(classification_report(y_test, y_pred1))
y_pred2 = KNclas_clf2.predict(X2_test)
print('\n---------- X2 -----------')
knn_eval1 = evaluate_scores(y_test,y_pred2)
print(classification_report(y_test, y_pred2))
y_pred3 = KNclas_clf3.predict(X3_test)
print('\n---------- X3 -----------')
knn_eval1 = evaluate_scores(y_test,y_pred3)
print(classification_report(y_test, y_pred3))
Again all three models have a worse f1 score than the dummy, the model trained on X2 is however close as it achieves the highest f1 score so far of 0.56. All three model outperformes the DummyClassifier on the overall accuracy and precision but the model trained on X2 again performes best. Any of these models are therefor more usable than the DummyClassifier. By further inspecting the precision for each category the precision of CLASE = 3 is relatively high for all models, but again the model trained on X2 does best. However this model also has a very low recall and f1 score for CLASE = 3. To further decide which model to choose we look at the confusion matrices.
# Confusion matrix
print('---------- X1 -----------')
confusion_ma(y_test,y_pred1)
print('---------- X2 -----------')
confusion_ma(y_test,y_pred2)
print('---------- X3 -----------')
confusion_ma(y_test,y_pred3)
By inspecting the confusion matrices we see that the model trained on X2 predicts the least amount amount of CLASE = 3 which is related to the low recall and thereby also f1 score for the CLASE. This model is therefor not choosen.
Of the two remaning models the KNN trained on X1 performes best. It by far predicts the highest amount of CLASE = 3 truly while only predicting 10% as CLASE = 1. The model also manages to predict CLASE = 1 better which can be seen on the f1 score and the confusion matrices. It is however out performed when looking at CLASE = 2, but, as described earlier, it is better to do worse in this class than the others. The best performing KNN is thus trained on the data set X1 where all relevant features are given with 1100 neighbors and uniform weights.
Random Forest
The random forrest classification is the last model which will be used for predicting. For the random forrest the following things are specified.
class_weight = 'balanced'- Adjusts class weights proportional to class frequencies.random_state = 42- To ensure same results every time the model is trained.n_estimators- For parameter for tuning. Is the number for trees in the forestmax_depth- For parameter tuning. Is the maximum depth of the tree.
A random forrest is fitted for X1_train.
from sklearn.ensemble import RandomForestClassifier
#Estimator
rfc = RandomForestClassifier(class_weight = 'balanced', random_state = 42)
#Parameters for the tuning
rfc_pars = {'n_estimators' : [9,11,13], 'max_depth' : [2,3,4]}
#GridSearchCV
rfc_clf1 = GridSearchCV(rfc, rfc_pars, cv=kf, scoring=scor)
rfc_clf1.fit(X1_train,np.ravel(np.ravel(y_train)))
#Results
print('Best parameters:',rfc_clf1.best_params_)
print('Best score:',rfc_clf1.best_score_)
A random forrest for X2_train is now fitted.
# #Parameters for the tuning
rfc_pars = {'n_estimators' : [8,9,10], 'max_depth' : [1,2,3]}
#GridSearchCV
rfc_clf2 = GridSearchCV(rfc, rfc_pars, cv=kf, scoring=scor)
rfc_clf2.fit(X2_train,np.ravel(np.ravel(y_train)))
#Results
print('Best parameters:',rfc_clf2.best_params_)
print('Best score:',rfc_clf2.best_score_)
Finally a random forrest for X3_train is fitted.
# #Parameters for the tuning
rfc_pars = {'n_estimators' : [13,14,15], 'max_depth' : [2,3,4]}
#GridSearchCV
rfc_clf3 = GridSearchCV(rfc, rfc_pars, cv=kf, scoring=scor)
rfc_clf3.fit(X3_train,np.ravel(np.ravel(y_train)))
#Results
print('Best parameters:',rfc_clf3.best_params_)
print('Best score:',rfc_clf3.best_score_)
We see that the models chooses different depths and amount of trees. With the fitted models we will look at the scores.
y_pred1 = rfc_clf1.predict(X1_test)
print('---------- X1 -----------')
rfc_eval1 = evaluate_scores(y_test,y_pred1)
print(classification_report(y_test, y_pred1))
y_pred2 = rfc_clf2.predict(X2_test)
print('\n---------- X2 -----------')
rfc_eval2 = evaluate_scores(y_test,y_pred2)
print(classification_report(y_test, y_pred2))
y_pred3 = rfc_clf3.predict(X3_test)
print('\n---------- X3 -----------')
rfc_eval3 = evaluate_scores(y_test,y_pred3)
print(classification_report(y_test, y_pred3))
The results reveal that all models do worse than the DummyClassifier on the overall f1 score, and also much worse compared to the logistic regression and KNN. However the overall accuracy and precision is much better than the DummyClassifier, and similar to the two other prediction methods. All three of these models are therefor better than the DummyClassifier. By inspecting the individual scores for each category all three random forrest models have similar f1 score for CLASE = 3 and that the recall equally is high. It is however noticable that the recall for CLASE = 2 is small for all three models. To decide which model performs best the confusion matrices are also inspected.
# Confusion matrix
print('---------- X1 -----------')
confusion_ma(y_test,y_pred1)
print('---------- X2 -----------')
confusion_ma(y_test,y_pred2)
print('---------- X3 -----------')
confusion_ma(y_test,y_pred3)
By inspecting the three confusion matrices it can be seen that the model trained on X3 manages to predict the largest share of CLASE = 3 correctly (86%). The model does however also predict 62% of CLASE = 2 and 33% of CLASE = 1 to be a CLASE = 3, this model is therefor not choosen.
Of the two remaning models the random forrest trained on X1 is choosen. This model is choosen as it only predicts 14% of CLASE = 1 to be CLASE = 3 where the model trained on X2 predict 20%. Furthermore the recall on CLASE = 2 is very low for the model trained on X2 which means that it only predicts 17% of CLASE = 2 correctly. Therefor the best random forrest is the model trained on X1 which is the data set with all relevant features. This model uses 11 trees and a maximum depth of 3.
With the best model for each of the classifiers selected the best model of these three is decided upon. To do this the scores and confusion matr
# Logistic regression
y_pred_logreg = logreg_clf3.predict(X3_test)
print('\n---------- Logistic regression, trained and tested on X3 -----------\n')
logreg_eval3 = evaluate_scores(y_test,y_pred_logreg)
print(classification_report(y_test, y_pred_logreg))
#confusion_ma(y_test,y_pred_logreg)
# KNN
y_pred_knn = KNclas_clf1.predict(X1_test)
print('\n---------- KNN, trained and tested on X1 -----------\n')
knn_eval1 = evaluate_scores(y_test,y_pred_knn)
print(classification_report(y_test, y_pred_knn))
#confusion_ma(y_test,y_pred_knn)
# Random forrest
y_pred_rfc = rfc_clf1.predict(X1_test)
print('\n---------- Random forrest, trained and tested on X1 -----------\n')
rfc_eval1 = evaluate_scores(y_test,y_pred_rfc)
print(classification_report(y_test, y_pred_rfc))
#confusion_ma(y_test,y_pred_rfc)
By comparing the scores for the three different classifiers we see that the logistic regression has the highest overall accuracy, precision and f1 score. The overall recall score is the second highest and is only 0.024 worse. For the individual scores for the different categories the f1 score for CLASE = 3 in the logistic regression is also best.
On the background of these observations the logistic regression which is trained on the dataset where all relevant features have been transformed using PCA (X3) is selected as the best predictor.
print('\n---------- Logistic regression, trained and tested on X3 -----------\n')
confusion_ma(y_test,y_pred_logreg)
print('\n---------- KNN, trained and tested on X1 -----------\n')
confusion_ma(y_test,y_pred_knn)
print('\n---------- Random forrest, trained and tested on X1 -----------\n')
confusion_ma(y_test,y_pred_rfc)
Prediction results¶
With a selected model the predicted results will now be plotted. Therefor a new dataframe containing the observed and predicted labels as well as the coordinates is defined. For the coordinates the earlier defined dataframe test_coor is used.
print(y_test.shape[0])
print(len(y_pred_logreg))
# Merge results and observed labels into a dataframe.
results = y_test
results['Predicted'] = y_pred_logreg
results = results.rename(columns={'CLASE':'True'})
# Merge the latitude and longitude
results['Longitude'] = test_coor['Longitude']
results['Latitude'] = test_coor['Latitude']
results.head()
The predicted labels are plotted together with the true labels.
colors = ['green', 'orange','dodgerblue']
fig, (ax1,ax2) = plt.subplots(1,2, figsize = (16,7))
# Scatter plot predicted
sns.scatterplot(x='Latitude', y='Longitude', hue='Predicted',
data=results, palette = colors, legend = 'full', ax = ax1)
ax1.set_title("Predicted labels")
# Scatter plot observed
sns.scatterplot(x='Latitude', y='Longitude', hue='True',
data=results, palette = colors, legend = 'full', ax = ax2)
ax2.set_title("True labels")
plt.show()
From the plots we see that there are too many households which are predicted as CLASE = 3 compared to the true labels and that most of these are gathered together in a cluster. This is as expected as we during the selection of the models put a high emphasis on that this CLASE was particularly important due to this being the most wealthy CLASE.
A figure showing the true and false labeled observations is also plotted.
results['Correctly_predicted'] = (results['Predicted'] == results['True'])
results.head()
colors = ['red','green']
plt.figure(figsize = (8,8))
sns.scatterplot(x='Latitude', y='Longitude', hue='Correctly_predicted',
data=results, palette = colors, legend = 'full')
plt.title("True and false predicted observations")
plt.show()
A area in the top of the figure is mostly dominated by true predicted observations. This is also the area as seen in the plot before where there was a large amount of true CLASE = 3.
fig,ax = plt.subplots(figsize=(15,5))
colors = ['red','green']
#ax = plt.subplot(1,2,2)
#sns.scatterplot(x="total_bill", y="tip", data=tips, ax= ax);
for i in range(1,4):
ax = plt.subplot(1,3,i)
sns.scatterplot(x='Latitude', y='Longitude', hue='Correctly_predicted',
data=results[results['True']==i], palette = colors, legend = 'full',ax=ax)
ax.set_title("Clase %i" %i)
The model is predicting CLASE 3 particularly good and does overall well. Clase 2 is not predicted very well. Maybe other household features than education would be better for prediciting.
However, to optimise the value for the real estate agency, the prediction should be presented by probability:
print(np.max(all_probs,0)) #Max probabilty
print(np.argmax(all_probs,0)) #Index for max probabilty
As an example the real estate agent sets real estate 11 for sale, the he gets the following decission support on who to target his advertisement:
print('This real estate has %.0f %% match with house hold type 1' %(logreg_clf3.predict_proba(X3_test[10:13])[0][0]*100))
print('This real estate has %.0f %% match with house hold type 2' %(logreg_clf3.predict_proba(X3_test[10:13])[0][1]*100))
print('This real estate has %.0f %% match with house hold type 3' %(logreg_clf3.predict_proba(X3_test[10:13])[0][2]*100))
In this case the real estate agent might benefit from mainly target advertising to type 1 but also target type 2.
Another example:
print('This real estate has %.0f %% match with house hold type 1' %(logreg_clf3.predict_proba(X3_test[1279:1280])[0][0]*100))
print('This real estate has %.0f %% match with house hold type 2' %(logreg_clf3.predict_proba(X3_test[1279:1280])[0][1]*100))
print('This real estate has %.0f %% match with house hold type 3' %(logreg_clf3.predict_proba(X3_test[1279:1280])[0][2]*100))
In this case the real estate agent can focus all the advertising efforts on type 3.
Conclusion¶
This report provides valuable insights for real estate agencies in Santiago de Chile. The models and insights derived in this report enable focused targeting of households based on eduction when advertising real estates. Predictions of the household type buying real estate in the city shows good results for predicting wealthy house holds with long education. These households are the ones that real estate agencies can earn most money on hence the most important. For the middle educated household the model is poor performing, potentially other household attributes than education would be better to classify here. For use in the business, the model can predict the probability of a real estate matching each of 3 house hold types. This allows the agencies to target relevant groups only. Additionally the market has been segmented with a fair match of household types. Together with a throughout exploratory analysis this report reveals key insights to understand the drivers for households to select real estate.
Contributions¶
Both members of the group have contributed equally and both members vouch for every section of the report. The construction of each section have been carried out as group work, both members have thereby contributed equally.